analytic  sciences  corp  reading  mass  f/g  8/6 

EFFICIENT  ESTIMATION  TECHNIQUES  FOR  INTEGRATED  GRAVITY  DATA  PRO— ETC(U) 
SEP  76  S W THOMAS#  W G HELLER  F19628-76-C-0014 

TASC-TR-680-1  AFGL-TR-76-0232  NL 


U.S.  DEPARTMENT  OF  COMMERCE 
National  Technical  Information  Service 


AD-A034  055 

EFFICIENT  ESTIMATION  TECHNIQUES  FOR 
INTEGRATED  GRAVITY  DATA  PROCESSING 

Analytic  Sciences  Corporation 
Reading/  Pennsylvania 

30  September  1976 


J 


AFGL-TR-76-0232 


EFFICIENT  ESTIMATION  TECHNIQUES  FOR  INTEGRATED  GRAVITY  DATA  PROCESSING 


Stephen  W.  Thomas 
Warren  G.  Heller 


The  Analytic  Sciences  Corporation 
Six  Jacob  Way 

Reading,  Massachusetts  01867 


Final  Report 

15  August  1975  - 31  July  1976 


Approved  for  Public  Release;  Distribution  Unlimited 


MR  FORCE  GEOPHYSICS  LABORATORY 

AIR  FORCE  SYSTEMS  COMMAND 

UNITED  STATES  AIR  FORCE 

Hanscom  AFB,  Massachusetts  01731 


REPRODUCED  BY 


NATIONAL  TECHNICAL 
INFORMATION  SERVICE 

U.S.  DEPARTMENT  OF  COMMERCE 
SPRINGFIELD,  VA.  221t>l 


'J2v 


4 


Qualified  requestors  may  obtain  additional  copies  from  the 
Defense  Documentation  Center.  All  others  should  apply  to 
the  National  Technical  Information  Service. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  of  this  PACE  flFhen  Doifi  Fm.i.d) 


\ REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

I.  REPORT  NUMBER 

AFGL-TR-76-0232 

2 GOVT  ACCESSION  NO. 

none 

J RECIPIENT'S  CATALOG  NUMBER 

|n  TITLE  (Mid  Submit) 

! Efficient  Estimation  Techniques  for 
Integrated  Gravity  Data  Processing 

s.  type  of  report  e perioo  covered 

Final 

8/15/75  - 7/31/76 

6 PERFORMING  ORG.  REPORT  NUMBER* 

TR-680-1 

|7  author^ 

6 CONTRACT  OR  GRANT  NUMBE Rf ») 

j Stephen  W.  Thomas 
1 Warren  G.  Heller 

F19628-76-C-0014 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

The  Analytic  Sciences  Corporation 
Six  Jacob  Way 

! Reading.  Massachusetts  01867 

10.  program  f.lement.  project,  task 

AREA  * WORK  UNIT  NUMBERS 

62101F 

76000301 

III.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Air  Force  Geophysics  Laboratory 
: Hanscom  AFB,  Massachusetts  1731 
\ Monitor/Geor^e  Hadgigeorge/LWG 

12  REPORT  DATE 

30  September  1976 

13  number  of  pages 

77 

i 14  MONITORING  AGENCY  NAME  5 ADDRESSfif  dlllermnt  from  Controlling  Office) 

is.  SECURITY  class.  (Of  thle  report) 

UNCLASSIFIED 

— - 

ISa  DC  CL  ASSl  T|  CATION  DOWNGRADING 
SCHEDULE 

16  DISTRIBUTION  statement  (of  thf.  R.port) 


! Approved  for  Public  Release  - Distribution  Unlimited 


17  DISTRIBUTION  STATEMENT  (of  the  abetrucl  enterod  In  Block  20.  If  different  from  Report ) 

ie  supplement  ary  notes 
i 

i 

) 19  KEY  WORDS  (Continue  on  reverre  si  dr  if  necessary  and  Identify  by  block  number) 

I Gravity  Disturbance,  Least-squares  Estimation,  Multisensor, 
Frequency  Domain,  Efficient  Computation 

| 20  ABSTRACT  (Continue  on  revereo  elde  If  neceeeary  nnd  Identify  by  block  number) 
I 

A mathematical  technique  is  developed  for  the  combined 

!use  of  various  types  of  gravimetric  data  in  the  estimation  of 
such  gravity  quantities  as  the  disturbance  potential  and  the 
gravity  disturbance  vector.  Modern  frequency  domain  techniques 
are  introduced  which  substantially  reduce  the  computer  process- 
ing requirements  associated  with  least-squares  gravity  quantity 
estimation.  This  approach  represents  an  extension  of  existing 


DD  ij?nM73  1 473  EDITION  OF  I NOV  *5  IS  OBSOLETE  UNCLASSIFIED 


I 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  Data  Lntemd) 


StCUNlTV  CLASSIFICATION  OF  THIS  NAOEf»Tl.n  Omit  Bnlm ,,d) 


J 


THE  ANALYTIC  SCIENCES  CORPORATION 


TABLE  GF  CONTENTS 

Page 
No . 


» List  of  Figures  iv 

List  of  Tables  v 

1.  INTRODUCTION  1-1 

1.1  Background  and  Purpose  1-1 

1.2  Scope  1-3 

1.3  Technical  Approach  1-5 

1.4  Overview  1-8 

2.  LEAST  SQUARES  ESTIMATION  OF  GRAVITY  QUANTITIES  2-1 

2.1  Statement  of  the  Problem  2-1 

2.1.1  Example  Problem  Formulation  2-2 

2.1.2  The  Role  of  Statistical  Models  2-4 

2.2  Least-Squares  Estimation  and  Error  Analysis  2-5 

2.3  Statistical  Gravity  Models  2-8 

3.  FREQUENCY  DOMAIN  ESTIMATION  WITH  SINGLE-SENSOR  DATA  3-1 

3.1  The  Frequency  Domain  Algorithm  3-1 

3.2  Relationship  With  Least-Squares  Estimation  3-4 

3.3  Numerical  Simulation  3-9 

3.4  Edge  Effect  Analysis  3-15 

4.  MULTISENSOR  FREQUENCY  DOMAIN  ESTIMATION  4-1 

4.1  Problem  Formulation  4-1 

4.2  Sequential,  Sensor-By-Sensor  Estimation  4-3 

4.3  Sensor-By-Sensor  Estimation  in  the  Frequency  Domain  4-6 

4.4  Multisensor  Numerical  Simulation  4-11 

5.  CONCLUSIONS  AND  RECOMMENDATIONS  5-1 

APPENDIX  A - ASYMPTOTIC  OPTIMALITY  OF  FINITE,  FREQUENCY 
k DOMAIN  ESTIMATION  A-l 

APPENDIX  B - MULTISENSOR,  DISCRETE  WIENER  FILTERING  B-l 

APPENDIX  C - DERIVATION  OF  FREQUENCY  DOMAIN  STATISTICS  USING 

THE  SAMPLING  THEOREM  C-l 

REFERENCES  R-l 


iii 


THE  ANALYTIC  SCIENCES  CORPORATION 


Figure 

No. 

1.1-1 

1.2-1 

1.3- 1 

1.3- 2 

2.3- 1 

3. 1- 1 

3.2- 1 

3.3- 1 

3.3- 2 

3.3- 3 

3.4- 1 

4.3- 1 

4.3- 2 

4.4- 1 


LIST  OF  FIGURES 


Integrated  Multisensor  Gravity  Data  Processing 

One-Dimensional,  Multisensor  Gravity  Survey 

Least  Squares  Estimation  and  Covariance  Analysis 

Relationships  Between  Direct  Least-Squares  Estimation 
and  the  Frequency  Domain  Approach 

Geometry  of  Round  Earth,  Attenuated  White  Noise  Model 
Single  Sensor  Frequency  Domain  Estimation 
The  Generalized  Wiener  Filter 

Projected  Gravity  Anomaly  Estimation  Error  for  an 
Airborne  Survey  With  Gradiometry  Only 

Power  Spectra  for  Airborne  Gradiometry 

Gravity- Anomaly/Gravity-Gradient  Spectral  Coherency 
Functions  for  Airborne  Gradiometry 

Estimation  Error  Sensitivity  to  Edge  Effects 

Sensor-By-Sensor  Frequency  Domain  Estimation 

Sampling  Theorem  for  Finite  Transforms 

Projected  Gravity  Anomaly  Estimation  Error  for  a 
Combined  Gradiometric/Gravimetric  Survey 


Page 

No. 

1-2 

1-4 

1-6 

1- 7 

2- 10 
3-3 
3-7 

3-12 

3-13 

3-14 

3- 17 

4- 8 
4-9 

4-13 


iv 


THE  ANALYTIC  SCIENCES  CORPORATION 


LIST  OF  TABLES 

Table 

No. 

2.2-1  Vector  Notation  for  Multisensor  Estimation 

3.1-1  Comparative  Computer  Processing  Times  for  Direct, 
Least  Squares  Estimation  and  Frequency 
Domain  Estimation 

4.4-1  Projected  Multisensor  Gravity  Survey  Accuracies 


Page 

No. 

2- 5 

3- 2 

4- 14 


THE  ANALYTIC  SCIENCES  CORPORATION 


1.  INTRODUCTION 

1.1  BACKGROUND  AND  PURPOSE 

Corresponding  to  the  need  for  more  accurate  naviga- 
tion and  guidance  systems,  there  has  been  an  increasing  in- 
terest in  detailed  knowledge  of  the  Earth's  gravity  field. 
This  interest  has  led  to  the  development  of  new,  sophisti- 
cated sensor  technology  to  augment  the  more  traditional 
gravimetric  and  astrogeodetic  techniques.  For  example, 
SEASAT-A  will  soon  be  in  orbit  providing  highly  accurate 
over-ocean  measurements  of  the  undulation  of  the  geoid. 
Similarly,  gravity  gradiometers  capable  of  providing  high 
quality  data  from  a moving-vehicle  survey  are  now  under 
development.  Satellite-to-satellite  tracking  techniques 
are  another  potential  data  source.  Thus,  large  amounts  of 
"multisensor  gravity  data"  will  be  available  in  the  near 
future . 


The  efficient  use  of  abundant  multisensor  gravity 
data  presents  new  problems  in  data  processing  and  analysis. 
Navigation  and  guidance  systems  respond  to  many  different 
gravity  disturbance  wavelengths,  depending  on  factors  such 
as  vehicle  speed  and  inertial  navigation  system  (Schuler 
frequency)  time  constants.  The  various  sensors  offer  differ- 
ing spectral  sensitivities  to  the  disturbances  in  the  Earth's 
gravity  field.  Gravity  gradient  data  best  describe  the  high 
frequency,  fine-grained  components  of  the  gravity  disturbance. 
Satellite  altimetry  and  sate) lite-to-satellize  tracking  tech- 
niques are  most  responsive  to  lower  frequency,  regional  gra- 
vity features.  Conventional  gravimetry  provides  still 
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another  data  source  with  an  intermediate  spectral  response. 

To  develop  a complete,  useful  description  of  the  Earth's 
gravity  field  requires  the  integrated  processing  of  data  from 
these  and  other  sensors  (Fig.  1.1-1).  Efficient  methods  for 
the  solution  of  this  multisensor,  gravity  data  processing 
problem  are  described  in  this  report. 


Mara 
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Figure  1.1-1  Integrated  Multisensor 
Gravity  Data  Processing 


This  report  extends  the  earlier  work  by  TASC  and 
others  on  gravity  data  processing  by  statistical  least- 
squares  techniques  (Refs.  1-5).  The  least-squares  approach 
has  certain  advantages.  Given  appropriate  statistical  grav- 
ity models,  this  methodology  provides  the  means  to  compute 
optimal  gravity  quantity  estimates  based  on  data  of  quite 
general  types.  The  accompanying  capability  for  statistical 
estimation  error  analysis  provides  the  method  with  a second 
attractive  feature.  Of  course,  the  problem  of  obtaining 
and  verifying  the  necessary  statistical  gravity  models  is  not 
a trivial  one.  While  not  a major  theme  in  this  report, 
research  in  this  important  problem  area  is  continuing  (Refs. 
6-11). 


Given  the  availability  of  appropriate  statistical 
models,  the  key  issue  in  least-squares  processing  of  multi- 
sensor gravity  data  is  the  development  of  techniques  to 
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mitigate  severe  computer  processing  requirements.  Even  for 
amounts  of  data  which  are  of  moderate  size  for  typical  grav- 
ity applications,  severe  computer  time  requirements  and  asso- 
ciated numerical  inaccuracies  make  the  method  extremely 
unattractive.  Computational  efficiency  must  be  an  inherent 
part  of  a practical  gravity  data  processing  algorithm.  This 
report  describes  a technique  which  provides  a substantial 
reduction  in  the  computational  cost  of  the  least-squares 
method;  as  applied  to  gravity  estimation. 

1 . 2 SCOPE 

This  report  is  concerned  with  the  development  of 
mathematical  techniques  for  the  estimation  of  certain  gravity 
quantities  from  sensors  which  are  currently  available,  or 
which  may  be  available,  in  the  near  future.  The  method  devel- 
oped here  are  appropriate  for  a broad  range  of  gravity  data 
sources,  such  as  gravity  gradiometry,  satellite  altimetry, 
gravimetry,  and  satellite-to-satellite  tracking.  The  esti- 
mation methods  of  this  report  are  constructed  under  the 
assumption  that  multisensor  data  is  available  on  one-dimen- 
sional tracks,  although  the  approach  used  here  is  applicable 
to  more  general  survey  geometries.  The  altitudes  of  the  data 
tracks  corresponding  to  different  sensors  may  differ,  as  may 
the  altitude  of  the  survey  track  along  which  estimates  are  to 
be  computed.  Numerical  results  are  specifically  presented 
for  gradiometry,  gravimetry,  and  satellite  altimetry.  Typical 
quantities  for  which  estimates  are  of  interest  are  gravity 
disturbance  potential  and  gravity  anomaly.  The  gravity  sur- 
vey configuration  treated  in  this  report  is  depicted  schema- 
tically in  Fig.  1.2-1.  The  problem  under  consideration  is  to 
obtain  efficient  gravity  quantity  estimates  based  on  moderate- 
to-large  amounts  of  multisensor  data. 
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A new  computationally-ef f icient  gravity  data  process- 
ing technique  is  described  which  employs  finite-dimensional 
frequency  domain  techniques.  This  approach  substantially  re- 
duces the  computer  processing  requirements  rroducing  mini- 

mum variance  gravity  quantity  estimates.  This  improvement 
in  efficiency  is  of  practical  importance  because  of  the  severe 
computational  requirements  of  existing  gravity  quantity  esti- 
mation (least-squares  collocation)  techniques.  A typical 
gravity  data  processing  problem  of  moderate  size  (10.000  data 
points)  might  require  several  days  of  computer  time  with  the 
existing  least-squares  collocation  methodology.  The  same 
problems  can  be  handled  in  less  than  a minute  with  frequency 
domain  methods. 

Numerical  examples  are  provided  which  demonstrate  the 
performance  of  the  frequency  domain  method  on  gravity  survey 
problems  of  the  type  described  above.  Used  in  conjunction 
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with  a new  TASC-developed  statistical  gravity  model  that  is 
particularly  attractive,  the  method  provides  a fast,  versa- 
tile gravity  quantity  estimation  tool. 


1.3  TECHNICAL  APPROACH 

The  technique  of  minimum  variance  estimation  provides 
the  general  framework  for  optimal  gravity  quantity  estimation 
as  discussed  in  this  report.  The  term  least-squares  colloca- 
tion is  often  applied  to  this  methodology  in  connection  with 
geophysical  applications.  This  technique  is  fundamentally 
statistical  in  nature,  and  relies  on  a statistical  model  in 
the  form  of  covariance  relationships  to  provide  the  connec- 
tion between  the  data  and  the  quantity  to  be  estimated,  and 
to  provide  a characterization  of  noise  in  the  data.  For 
example,  if  airborne  gravity  gradiometry  measurements  are 
used  to  estimate  the  gravity  anomaly  at  the  earth's  surface, 
a covariance  matrix  must  implicitly  characterize  the  physi- 
cal upward  continuation  of  the  gravity  disturbance  to  the 
aircraft  altitude.  Other  covariance  matrices  must  character- 
ize the  quality  of  any  prior  knowledge  of  gravity  anomaly  in 
the  survey  region,  as  well  as  the  random  noise  in  the  gradiom- 
eter  itself.  The  least-squares  estimation  equations  specify 
the  operations  that  must  be  carried  out  with  these  covariance 
matrices  to  form  optimal  (minimum  mean-squared  error)  gravity 
quantity  estimates  based  on  the  data.  An  error  covariance 
analysis  based  on  these  equations,  and  containing  matrix 
operations  of  similar  form,  predicts  the  squared  error  in 
the  resulting  estimates.  These  relationships  are  diagrammed 
in  Fig.  1.3-1. 

Frequency  domain  methods  provide  a substantial 
reduction  in  the  numerical  cost  of  performing  the  matrix 
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Figure  1.3-1 


Least  Squares  Estimation  and 
Covariance  Analysis 


operations  inherent  in  the  least  squares  estimation  equations 
(Refs.  12-15).  The  key  to  this  result  is  that  covariance 
matrices  of  a certain  general  structure  become  diagonal  in 
the  frequency  domain.  The  crucial  operations  of  matrix  in- 
version and  matrix  multiplication  then  require  substantially 
fewer  computer  operations  (multiplications,  additions,  retrie- 
vals from  storage)  than  in  the  full-matrix  case.  This  reduc- 
tion in  computational  cost  more  than  compensates  for  the  cost 
of  transformations  to  and  from  the  frequency  domain.  The  con- 
ceptual relationship  between  the  direct  application  of  the 
least-squares  estimation  equations  and  the  frequency  domain 
approach  is  illustrated  in  Fig.  1.3-2. 


The  requirements  for  the  applicability  of  the  fre- 
quency domain  method  are  simply  related  to  the  mathematical 
structure  of  the  covariance  matrices  involved.  In  terms  of 
the  gravity  quantity  estimation  problem  itself,  it  is  re- 
quired that  the  data  be  available  at  uniformly  spaced  points, 
although  interpolated  data  is  acceptable,  and  that  the  survey 
track  be  long  compared  to  the  distance  over  which  the  rele- 
vant gravity  quantity  remains  highly  correlated.  These 
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Squares  Estimation  and  the  Frequency 
Domain  Approach 

requirements  seem  reasonable  for  the  gravity  survey  types  con- 
sidered in  this  report. 

0 

The  extension  of  the  techniques  of  this  report  to 
gravity  quantity  estimation  on  a two-dimensional  survey 
region  seems  appropriate  under  requirements  similar  to  those 
above.  Here  the  need  for  uniformly  spaced  data  translates 
into  a requirement  for  uniform  data  intervals,  measured  as 
spherical  angles,  along  both  axes  of  a two-dimensional  grid. 
Interpolated  or  averaged  data  remains  a possibility  for  the 
two-dimensional  problem.  Under  these  conditions,  an  appli- 
cation of  the  two-dimensional,  finite  Fourier  transform  would 
lead  to  a diagonalization  of  the  estimation  equations,  analog- 
ous to  the  one-dimensional  case.  Extension  of  th  frequency 
domain  methods  of  this  report  offers  all  of  the  computational 
advantages  associated  with  the  one-dimensional  procedure. 


GRAVITY  QUANTITY 
ESTIMATE.  £ 


INVERSE  TRANSFORM 
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1.4  OVERVIEW 

The  mathematical  statement  of  the  minimum  variance 
gravity  quantity  estimation  problem  and  its  general  solution 
are  presented  in  Chapter  2 . The  application  of  these  general 
expressions  to  a given  problem  requires  specific  statistical 
models,  and  it  is  shown  how  these  may  be  obtained  from  a new 
TASC-developed  gravity  disturbance  potential  model. 

Chapter  3 outlines  the  derivation  of  the  finite 
frequency  domain  estimation  method  for  single-sensor  data. 

A covariance  simulation  demonstrates  the  application  of  the 
method  to  estimation  of  the  gravity  anomaly  using  airborne 
gravity  gradiometer  data.  Further  numerical  work  evaluates 
the  impact  of  finite  track  length  approximations  on  the  util- 
ity of  the  frequency  domain  method. 

Frequency  domain  methods  for  least-squares  estima- 
tion with  multisensor  data  are  considered  in  Chapter  4 . In 
some  important  situations  with  multisensor  data  the  frequency 
domain  equations  are  no  longer  diagonal,  but  retain  useful 
structure.  A sequential,  one-sensor-at-a-time  method  is 
presented  which  recovers  the  diagonal  structure  in  these 
cases.  The  multisensor  frequency  domain  technique  is  demon- 
strated with  error  covariance  analyses  for  surveys  using 
airborne  gradiometry,  satellite  altimetry,  and  surface  gravi- 
metry in  combination  for  the  estimaton  of  the  gravity  anomaly. 

Chapter  5 reviews  the  results  of  this  report  and 
discusses  further  applications  and  extensions  of  frequency 
domain  techniques  in  multisensor  gravity  data  processing. 
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2.  LEAST  SQUARES  ESTIMATION  OF  GRAVITY  QUANTITIES 

2.1  STATEMENT  OF  THE  PROBLEM 

12  N 

Measurements,  z , z , ...,  z of  variables  correspond- 
ing to  N different  sensors,  which  may  be,  for  example,  gravity 
gradiometry,  satellite  altimetry,  or  conventional  gravimetry 
are  available  along  one  dimensional  data  grids.  These  observa- 
tions contain  additive  random  instrumentation  noise.  The  prob- 
lem is  to  compute  optimal  estimates  of  a gravity  disturbance 
quantity  x,  which  may  typically  be  the  gravity  anomaly,  verti- 
cal deflections,  or  disturbance  potential,  at  grid  points  on 
or  above  a one-dimensional  survey  track  (Fig.  1.2-1).  Note 
that  the  data  and  survey  tracks  need  not  all  be  at  one  alti- 
tude, so  that  essentially  a two-dimensional  survey  geometry 
is  being  considered.  The  estimates  for  x are  to  be  optimal 
with  respect  to  a statistical  model  which: 

• Is  consistent  with  the  physical  constraints 
between  x and  z1 , ....  zN  (e.g.,  upward/ 
downward  continuation). 

• Contains  prior  statistical  knowledge  or 
assumptions  concerning  x,  z1,  ...,  zN. 

• Describes  the  characteristics  of  the  measure- 
ment noise  present  in  the  data. 

The  estimate  x sought  for  x is  to  be  optimal  in  the  sense  of 
minimizing  the  mean  squared  error  Et(x-x)2],  where  E [ • ] denotes 
the  ensemble  expectation  operator. 

The  point  sensor  "measurements"  may  themselves  be 
the  result  of  the  preprocessing  of  some  raw  data.  In  the 
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case  of  airborne  gradiometry,  for  example,  the  data  might  be 
obtained  as  a sampled  moving  window  average  of  a continuous 
instrument  measurement.  Surface  gravimetry  data  might  be  the 
result  of  interpolating  or  averaging  data  from  some  more 
complicated  survey  geometry. 

The  relationships  among  the  mathematical  estimation 
problem,  the  underlying  physical  variables,  and  the  statis- 
tical models  involved  are  best  described  with  a concrete 
example.  The  example  presented  here  will  be  followed  up  in 
subsequent  sections,  culminating  with  a numerical  example  in 
Section  3.3. 

2.1.1  Example  Problem  Formulation 


It  is  assumed  that  a single  element  of  the  gravity 

gradient  tensor  is  to  be  measured  via  airborne  gradiometry, 

and  that  the  information  gathered  is  to  be  used  to  estimate 

the  gravity  anomaly  along  the  aircraft's  flight  path.  The 

gravity  gradient  at  the  aircraft  altitude  is  measured  at 

points  along  a regular  grid.  The  gravity  anomaly  is  to  be 

estimated  at  the  points  of  a similar  grid  on  the  surface. 

These  quantities  are  most  easily  handled  in  vector  form. 

The  gridded  values  of  the  gravity  anomaly  and  the  observed 

+ 

gradient  are  denoted 


x = (x-p  . . . ,xn) 


z = (z2 zn) 


(gravity  anomalies) 


(observed  gradients) 


(2.1-1) 


respectively.  For  simplicity,  it  has  been  assumed  that  the 
two  sampling  grids  are  of  the  same  lengths  and  spacing,  so 
that  x and  z are  of  the  same  dimension. 

The  superscript  ( ) denotes  the  transpose  operation. 


2-2 


THE  ANALYTIC  SCIENCES  CORPORATION 


The  required  statistical  models  for  x and  z may 
be  divided  into  a statistical  gravity  model,  which  describes 
the  gravity  anomaly  and  its  relationship  to  the  true,  or  noise- 
free,  value  of  the  anomalous  gradient,  and  a model  for  measure- 
ment noise  in  the  observed  values  of  the  anomalous  gradient. 

It  is  convenient  to  express  these  observations  as 

z = £ + v ( 2 . 1-2 ) 

where  £ contains  the  values  of  the  "true"  anomalous  gradient, 
and  v consists  of  random  measurement  noise. 

The  statistical  gravity  model  consists  of  the  dis- 
crete auto-  and  cross  covariance  functions 

»xx<k>  = E(xj  xj+k> 

*yy<k)'*  E(yj  yj+k> 

*xy(k)  - E(xj  yj+k)  (2.1-3) 

for  1 < j,(j+k)  1 n.  The  measurement  noise  v is  modeled  as 
being  statistically  uncorrelated  with  x and  £,  with  auto- 
correlation given  by 

*vv(k)  ■=  E(v.  vj+k)  (2.1-4) 

All  three  variables  x,  £,  z have  zero  mean. 

The  specification  of  the  statistical  models  completes 
the  definition  of  the  mathematical  estimation  problem  for  the 
example  at  hand.  The  specifics  of  the  computation  of  opti- 
mal estimates  will  be  considered  in  Section  2.2. 
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2.1.2  The  Role  of  Statistical  Models 

Apart  from  the  mechanics  of  computing  gravity  anomaly 
estimates,  the  role  of  statistical  models  in  the  example  just 
presented  merits  discussion  at  this  point . Perhaps  the  most 
important  factor  in  a gravity  quantity  estimation  problem  of 
this  type  is  the  set  of  physical  consistency  relationships 
that  link  the  observable  variables  with  those  to  be  estimated. 
These  relationships  will  be  properly  accounted  for  in  the  esti- 
mated process  if,  and  only  if,  they  are  embodied  in  the  asso- 
ciated statistical  gravity  model.  For  the  airborne  gradiometry 
example  this  means  that  the  operations  of  upward  continuation 
and  differentiation  which  relate  the  gravity  anomaly  and  the 
observed  anomalous  gradient  must  be  implicit  in  Eq.  2.1-2.  For 
this  reason,  the  idea  of  self-consistency  in  the  construction 
of  gravity  models,  as  considered  in  Refs.  8 and  9,  is  very  im- 
portant to  the  estimation  process. 

The  estimation  problem  described  in  this  section  c?„n 
be  extended  by  introducing  uncertain,  but  nonrandom,  parame- 
ters into  the  measurement  expression,  Eq.  2.1-2: 

z = + v + y.  (2.1-5) 

where  S is  a vector  of  unknown  parameters  and  H is  a matrix 
with  compatible  dimensions.  The  term  Hj)  represents  "system- 
atic effects"  present  in  the  data.  This  can  be  useful  in 
modeling  effects  such  as  sensor  biases  or  drifts  which  lead 
to  trends  in  the  data.  The  estimation  of  £ and  £ from  mea- 
surements of  the  form  Eq.  2.1-5  has  been  considered  using  the 
method  of  least-squares  collocation  (Refs.  3,  4).  This  treat- 
ment leads  to  estimation  expressions  similar  to  those  pre- 
sented in  Section  2.2,  which  are  based  on  Eq . 2.1-2.  The 
approach  to  least-squares  estimation  based  on  Eq . 2.1-2  is 
used  in  this  report  because  it  leads  more  directly  to  the 
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frequency  domain  methodology.  The  explicit  application  of 
frequency  domain  methods  to  models  using  Eq.  2.1-5  would  be 
a useful  continuation  of  the  present  work. 


2.2  LEAST-SQUARES  ESTIMATION  AND  ERROR  ANALYSIS 

The  solution  to  the  least-squares  estimation  prob- 
lem is  most  easily  presented  in  vector-matrix  notation.  The 
column  vector  notation  used  in  the  sequel  is  defined  in 
Table  2.2-1.  The  ordering  of  multisensor  data  within  £ and 
z is  left  unspecified  because  it  does  not  have  impact  on 
the  general  estimation  expressions  presented  in  this  section. 


TABLE  2.2-1 

VECTOR  NOTATION  FOR  MULTISENSOR  ESTIMATION 


VECTOR 


DEFINITION 

DIMENSION 

Gridded  values  of  quan- 
tity to  be  estimated 

n 

True  values  of  observed 
variables  (all  sensors) 
at  grid  points 

m 

Random  measurement  noise 

m 

Observed  data  (z  = £ + v) 

m 

The  information  provided  by  the  statistical  models 
associated  with  the  above  variables,  typically  available  in 
the  form  of  Eq.  2.1-3,  2.1-4  is  conveniently  expressed  in 
covariance  matrix  form.  The  p*q  matrix 


cab  - E<i  S > 
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denotes  the  covariance  matrix  of  random  (zero  mean)  vectors 
a and  b of  dimension  p and  q,  respectively.  The  statistical 
covariance  model  associated  with  the  estimation  problem  con- 
sists of  the  matrices  C , C , C , C , which,  together  with 

aa  yy  xy  vv 

the  assumption  that  £ and  v are  independent,  yields  the  data 
covariance  matrix 


C + C 
yy  vv 


(2.2-2) 


The  expressions  for  minimum  variance  estimation  and 
error  covariance  analysis  are  well  known  (Refs.  16,17).  The 
optimal  estimate  for  x,  based  on  the  model  provided  by  Eqs. 
2.1-2  through  2.1-4,  is 


~ -1 
x = C C z 
— xz  zz  — 


(2.2-3) 


This  expression  utilizes  the  assumption  that  x and  z have 
zero  a priori  means.  If  this  is  not  the  case,  the  vectors 

A 

x and  z in  Eq . 2.2-3  must  be  replaced  by  their  respective 
departures  from  the  means  of  x and  z.  A statistical  pre- 

A 

diction  of  the  estimation  error  associated  with  x as  computed 
by  Eq.  2.2-3  may  be  obtained  as  the  covariance  matrix  of 
the  estimation  error,  e = x - x.  The  latter  error  matrix 
may  be  computed  from 


C - C C-1  C 
xx  xz  zz  zx 


(2.2-4) 


Although  the  matrix  C gives  a complete  statistical 
picture  of  the  errors  in  x and  their  correlations,  it  is  use- 
ful to  relate  the  overall  estimation  error  in  a given  situa- 
tion to  a single  quantity.  The  rms  estimation  error,  o , 
defined  by 
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e 


1 

n 


E e[" 
:=1  L 


(xkA 


(2.2-5) 


will  be  used  extensively  in  the  sequel.  For  the  one- 

dimensional  gravity  survey  grid  considered  in  this  report, 

2 

og  may  be  interpreted  as  the  average  squared  estimation 
error  along  the  survey  track.  It  is  readily  calculated  from 

Cee  using* 


a2  = - tr(C  ) 
e n eey 


(2.2-6) 


While  Eqs.  2.2-3,  2.2-4  provide  a very  general 

vehicle  for  the  calculation  of  optimal  gravity  estimates, 

in  problems  of  even  moderate  size  they  present  computer 

processing  problems  that  cannot  be  overlooked.  Both 

expressions  contain  C_1 , the  inverse  of  a matrix  whose 

zz 

dimension  equals  the  number  of  data  elements,  m.  The  most 

0 

efficient  general  algorithms  for  producing  the  required 
inverse  (which,  for  efficiency,  would  actually  involve 
solving  a system  of  linear  equations)  call  for  a number  of 

3 

numerical  operations  proportional  to  m . This  has  a very 
strong  limiting  effect  on  the  amount  of  data  for  which  it 
is  economical  to  apply  Eqs.  2.2-3,  and  2.2-4.  For  example, 
on  a modern,  fast  computer  15  minutes  would  be  sufficient 
to  calculate  the  vector  C_1z  for  m=1000,  but  approximately 

ZZ — 

10  days  would  be  required  for  m=10,000.  Based  as  much  on 
numerical  accuracy  problems  as  on  consumption  of  computer 
time,  a figure  of  1000  data  elements  (m=1000)  may  be 
taken  as  a practical  upper  bound  for  the  successful 
direct  application  of  Eqs  2.2-3  and  2.2-4. 


\ 


♦The  symbol  tr(C  ) denotes  the  trace  of  the  matrix  C 

ee  ee 

i.e.,  the  sum  of  the  diagonal  elements  of  C 

ee 


l 
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The  airborne  gradiometry  example  of  the  last  section 
helps  to  put  the  numerical  processing  problem  into  perspective 
with  the  physical  availability  of  data.  If  gravity  gradient 
data  were  recorded  at  1 nm  intervals,  then  a survey  flight  of 
moderate  length  might  yield  data  at  1000  locations  along  the 
survey  track.  If  only  one  element  of  the  gravity  gradient 
tensor  is  measured,  then,  by  taking  great  care  with  numerical 
algorithms,  Eqs.  2.2-3  and  2.2-4  might  be  used  directly.  If, 
however,  more  than  one  of  the  five  independent  gradient  ele- 
ments is  measured  at  each  location,  then  Eqs.  2.2-3  and  2.2-4 
as  they  stand  are  probably  not  satisfactory  because  of  the 
computer  time  requirement  alone.  This  problem  is  further 
compounded  by  the  addition  of  data  from  any  other  sensors,  or 
of  data  from  other,  parallel  survey  tracks.  The  need  to  effi- 
ciently apply  Eqs.  2.2-3  and  2.2-4  to  gravity  estimation  prob- 
lems of  this  type  provides  the  motivation  for  the  frequency 
domain  approach  to  computation  developed  in  the  sequel. 


2.3  STATISTICAL  GRAVITY  MODELS 

The  optimal  estimation  expressions  require  statis- 
tical covariance  functions  for  the  gravity  disturbance  quan- 
tities involved.  As  discussed  earlier  in  the  chapter,  the 
choice  of  the  covariance  functions  used  for  this  purpose  must 
be  restricted  to  those  which  obey  the  appropriate  physical 
consistency  relations.  The  construction  of  gravity  disturb- 
ance models  has  been  considered  by  a number  of  authors 
(Refs.  7 to  11). 

Of  the  available  self-consistent  models,  those  of 
Refs.  9,  10,  and  11  are  the  only  ones  that  provide  covari- 
ance expressions  for  points  at  different  altitudes  without 
requiring  analytically  intractable  or  numerically  expensive 
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upward  continuation  operations.  The  models  of  Heller  (Ref.  11) 
and  Bellaire  (Ref.  10)  have  the  practical  advantage  that  they 
are  simple,  and  therefore  quite  tractable.  The  model  used  for 
the  purposes  of  this  report  is  the  asymptotic  form  of  the 
spherical  earth  "attenuated  white  noise"  (AWN)  model  developed 
in  Ref.  11.  This  model  is  mathematically  equivalent  to  the 
one  derived  in  Ref.  10  based  on  "flat  earth"  geometry.  As 
discussed  below,  this  asymptotic  form  of  the  AWN  model  is 
somewhat  more  generally  applicable  than  the  "flat  earth"  de- 
rivation of  Ref.  10  would  suggest,  and  yet  is  simpler  than 
the  more  general  spherical  earth  AWN  model. 

The  geometry  of  the  spherical  earth  AWN  model  is 
illustrated  in  Fig.  2.3-1.  Suppose  the  disturbance  poten- 
tial on  a spherical  shell  at  depth  D is  taken  to  be  uncorre- 
lated from  point  to  point,  i.e.,  white  noise.  Then  the 
appropriate  upward  continuation  integral  may  be  solved  analy- 
tically for  the  disturbance  potential  autocovariance,  4>TT(  ^ , 
rl’  r2^  between  points  at  radii  r^  and  r^  separated  by  a spher- 
ical angle  i p.  A useful  approximation  to  this  autocovariance 
function  may  be  obtained  by  taking  the  limit  as  the  ratio 
of  shell  depth  to  the  Earth's  radius,  and  shift  angle  approach 
zero.  The  resulting  asymptotic  form  of  the  AWN  model  distur- 
bance potential  autocovariance  function  is 


4>tt^  ^ 


4D2(2D+h1+h2)a£ 


~(2D+h1+h2)2  + t2J 


2-1372 


(2.3-1) 


where  t is  the  shift  distance  between  points  and  h^ , h2  are 
heights  above  the  earth's  surface.  A disturbance  potential 
model  of  more  general  form  can  be  created  by  superposing 
versions  of  Eq . 2.3-1  with  differing  parameters  (oT  and  D). 
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R171S4 


CORRELATION 

POINTS 


figure  2.3-1  Geometry  of  Round  Earth, 

Attenuated  White  Noise  Model 


The  result  of  such  a superposition  can  be  viewed  as  a model 
based  on  several  white  noise  shells  with  differing  intensities 
and  depths.  The  ability  of  this  procedure  to  represent  real 
data  is  indicated  in  Ref.  11,  where  a "three  shell"  model  has 
been  used  to  approximate  a worldwide  gravity  anomaly  auto- 
covariance function  with  less  than  1%  rms  error. 

It  is  important  to  observe  that  the  application  of 
Eq . 2.3-1  is  not  restricted  to  situations  where  flat  earth 
geometry  may  be  expected  to  be  a good  physical  approximation. 
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Analysis  presented  in  Ref.  11  shows  that  Eq.  2.3-1  is  a good 
approximation  to  its  spherical  analog  for  shell  depths,  D, 
less  than  several  hundred  miles,  with  no  restrictions  on 
altitude  above  the  surface.  The  reader  is  referred  to  Ref.  11 
for  a derivation  and  analysis  of  the  asymptotic  form  of  the 
AWN  model,  as  well  as  for  a demonstration  of  its  capability 
to  fit  real  data. 


Given  an  expression  for  the  disturbance  potential 
autocovariance  function,  the  required  covariance  function 
models  for  other  anomalous  gravity  quantities  may  be  obtained 
through  simple,  mathematical  operations.  The  undulation  of 
the  geoid  is  of  interest  for  gravity  survey  applications  in- 
volving satellite  altimetry.  The  necessary  statistical  model 
for  the  undulation,  N,  may  be  obtained  directly  from  Eq . 2.3-1 
using  Brun's  formula  (Ref.  18): 


N = |-  (2.3-2) 

go 

where  gQ  is  the  nominal  local  magnitude  of  the  gravity  vec- 
tor. Covariance  models  for  gravity  anomaly,  the  elements  of 
the  gravity  disturbance  vector,  and  the  elements  of  the  gravity 
gradient  tensor  are  required  in  survey  applications  using 
gravimetry  or  gravity  gradiometry.  These  quantities  are 
physically  related  to  the  disturbance  potential  through  the 
appropriate  derivative  operations  (Ref.  18).  It  follows  from 
the  elementary  properties  of  covariance  functions  (Ref.  19) 
that  self-consistent  models  for  these  quantities  can  be  pro- 
duced by  applying  corresponding  derivative  operations  to  the 
disturbance  potential  autocovariance  function.  These  are  tab- 
ulated in  Appendix  F of  Ref.  11.  Expressions  for  these  deri- 
vatives can  be  generated  by  an  automated  procedure,  such  as 
the  IBM  FORMAC  routine. 
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The  attenuated  white  noise  model  of  Eq . 2.3-1, 
together  with  Eq . 2.3-2  and  the  capability  for  automated 
differentiation,  provides  a readily  accessible,  self- 
consistent  statistical  gravity  model  of  great  flexibility 
All  of  the  statistical  models  for  the  numerical  examples 
in  this  report  were  generated  in  this  way. 


2-12 


THE  ANALYTIC  SCIENCES  CORPORATION 


3.  FREQUENCY  DOMAIN  ESTIMATION  WITH 

SINGLE-SENSOR  DATA 

3.1  THE  FREQUENCY  DOMAIN  ALGORITHM 

Frequency  domain  estimation  is  a fast  computational 
technique  for  carrying  out  least-squares  estimation  and  esti- 
mation error  analysis.  For  a broad  class  of  estimation  prob- 
lems it  produces  essentially  the  same  result  as  the  direct 
least-squares  procedure  of  Section  2.2,  but  at  a small  frac- 
tion of  the  computational  cost.  This  is  possible  because 
the  transformation  to  the  frequency  domain  allows  the  struc- 
tural regularities  already  present  in  the  estimation  expres- 
sions to  be  put  to  good  use.  For  the  simplest  case,  the 
least-squares  equations,  which  are  full  but  highly  regular 
in  their  usual  representation,  become  diagonal  in  their 
frequency  domain  form.  The  reduction  in  cost  of  carrying 
out  the  simplified  equations  more  than  compensates  for  the 
required  transformations  to  and  from  the  frequency  domain. 
Comparative  estimates  of  numerical  processing  times  (IBM 
370/165)  for  frequency  domain  and  direct  estimation  are  pre- 
sented in  Table  3.1-1.  When  the  number  of  data  points,  n, 
is  large,  the  time  required  for  the  direct  method  is  larger 
by  a factor  proportional  to  n /log  n. 

The  statistical  theory  of  spectral  representation 
(Refs.  20  - 22)  provides  elegant  motivation  for  the  use  of 
frequency  domain  estimation  methods.  This  body  of  mathemati- 
cal theory  is  concerned  with  the  statistical  characteristics 
of  Fourier-transformed  random  processes.  Let  x be  a random 

n-vector,  and  let  X denote  its  discrete  Fourier  transform 

* 

(DFT),  defined  as  the  column  vector  with  the  elements 
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TABLE  3.1-1 

COMPARATIVE  COMPUTER  PROCESSING  TIMES  FOR  DIRECT,  LEAST- 
SQUARES  ESTIMATION  AND  FREQUENCY  DOMAIN  ESTIMATION 


Number  of  Data  Points 

1000 

10,000 

50,000 

Direct  Estimation* 

14  min 

10  days 

~3  years 

Frequency  Domain 
Estimation 

<<  1 sec 

1.4  sec 

10  sec 

X(U)  = Xk+1  exp(2lis*)  (3.1-1) 

/n  k=0  K 1 V n / 

for  u)  = 0,...,n-l.  Suppose  that  x is  statistically  stationary 
in  the  sense  that 

E(xk  ’W  = ^xx(j)  : 1 1 j-(J+k)  1 n 

An  important  result  of  spectral  representation  theory  (Ref. 

21)  states  that  the  frequency  domain  elements  X(uj)  and  X(p) 
are,  within  an  "edge  effect"  approximation,  statistically 
uncorrelated  for  u>  f p . A similar  result  for  any  pair  of  ran- 
dom vectors,  x and  z,  provides  for  the  lack  of  correlation 
between  the  transforms  X(co)  and  Z(p)  with  gj  f p.  These  are 
the  key  results  which  motivate  frequency  domain  estimation 
techniques . 


The  application  of  spectral  representation  theory 
to  the  least-squares  estimation  problem  leads  to  an  algorithm 
which  is  a variation  on  the  classic  Wiener-Kolmogorov  filter 
(Ref.  23).  Let  x and  z be  defined,  following  the  notation 
of  Section  2.2,  as  the  vector  to  be  estimated  and  the  data 

♦ , 

The  letter  i denotes  the  complex  number  /-l. 

tProcessing  times  based  on  use  of  a direct  matrix  inversion 

procedure  to  solve  Eq.  2.2-3. 
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vector,  respectively.  The  spectral  representations  of  the 
elements  of  these  vectors,  X(w)  and  Z(w),  can  be  computed 
using  Eq.  3.1-1.  Under  regularity  conditions  corresponding 
to  Eq.  3.1-2,  the  values  of  these  transforms  at  any  one  freq- 
uency are  uncorrelated  with  their  values  at  any  other  freq- 
uency. This  implies  that  least-squares  estimation  of  X(u>) 
can  be  carried  out  based  only  on  the  data  at  the  same  freq- 
uency, Z(u>).  The  application  of  the  least-squares  equations, 
Eq.  2.2-3  and  2.2-4,  one  frequency  at  a time  leads  to 


XO) 


Z(w) 


(3.1-3) 


$ ((jj)  = <J>  (oj)- 

eev  ; xx 


*xz(a))  *zx(co) 
*zz(li)) 


(3.1-4) 


The  quantities  * (a >),  $ ( to ) , ancl  4>  (u>)  correspond  to  the 

XX  zz  xz 

appropriate  frequency  domain  variances,  and  covariances,  of 
X(io)  and  Z(uj).  The  expressions  in  Eq . 3.1-3  and  3.1-4  form 
the  basis  of  the  single-sensor,  frequency  domain  estimation 
algorithm.  The  overall  procedure  is  presented  in  Fig.  3.1-1 


ESTIMATION  ALGORITHM 


R441M 


ESTIMATION  ERROR  ANALYSIS 


CrtCMC*X 


Figure  3.1-1  Single  Sensor  Frequency  Domain  Estimation 
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The  requirements  for  the  applicability  of  the  fre- 
quency domain  estimation  procedures  are  based  primarily  on 
the  stationarity  condition,  Eq.  3.1-2.  The  statistical  gravity 
model  of  Section  2.3  is  of  the  stationary  type,  as  are  all 
of  the  models  proposed  in  Refs.  6-11.  For  x and  z based 
on  grids  with  uniform  spacing,  and  when  only  single  sensor 
data  is  used,  Eq.  3.1-2  will  be  satisfied  with  these  models. 

A more  complete  set  of  necessary  conditions  is  presented  in 
Section  3.2. 


Computation  with  the  frequency  domain  estimation 
method  is  facilitated  by  some  useful  numerical  devices.  The 
fast  Fourier  transform  (FFT)  algorithm  (Refs.  24,25)  can  be 
used  to  efficiently  carry  out  the  required  transforms  and 
inverse  transforms.  The  use  of  the  FFT  accounts  for  a part 
of  the  superior  speed  of  the  method.  The  statistics  <t>  ( co ) , 

XA 

$ (u>)>  and  4>  (u>)  are  computed  as  power  spectral  densities 

zz  xz 

based  on  the  matrices  C . C.  and  C „ which  define  the  ori- 

XX  zz  xz 

ginal  estimation  problem.  The  mathematical  definition  for 
these  power  spectra  is  provided  in  Appendix  A (Eq.  A-5). 
Finally,  the  rms  estimation  error,  a , can  be  computed  without 
performing  the  inverse  transform  of  Xee(u>)  by  using  a form  of 
Parseval's  identity  (Ref.  25) 


I *ee(<*>) 

u)-0  e 


(3.1-5) 


3.2  RELATIONSHIP  WITH  LEAST-SQUARES  ESTIMATION 

Generalized  Wiener  filter  theory  (Refs.  26,  27) 
provides  a framework  for  establishing  the  relationship  between 
direct  least-squares  estimation  and  the  frequency  domain  esti- 
mation method.  This  approach  reveals  the  frequency  domain 
method  simply  as  an  algebraically  transformed  version  of  the 


I 


B 
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standard  least-squares  procedure,  modified  for  computational 
efficiency.  One  result  of  this  treatment  is  insight  into  the 
asymptotic  equivalence  of  the  estimates  produced  by  the  direct 
least-squares  and  frequency-domain  methods.  Edge  effect 
approximations  cause  the  estimates  computed  with  the  frequency 
domain  method  to  depart  somewhat  from  the  exact,  least-squares 
solution.  Under  the  appropriate  conditions  the  departure  from 
the  least-squares  estimates  becomes  arbitrarily  small  as  sur- 
vey track  length  is  increased.  Another  important  application 
of  the  methodology  outlined  in  this  section  is  the  development 
of  frequency  domain  techniques  appropriate  for  multisensor 
problems,  as  discussed  in  Chapter  4. 

A block  diagram  of  the  generalized  Wiener  filter 
process  is  illustrated  in  Fig.  3.2-1.  This  technique  is 
based  upon  the  fact  that  least-squares  estimation  may  be 
equivalently  carried  out  under  any  appropriate  transformation. 
Such  a transformation  is  defined  by  applying  a nonsingular 
matrix  F via 


Z = Fz  , X = Fx  (3.2-1) 

and  computing  the  least-squares  estimate  of  X based  on  Z. 

The  matrix  F may  be  complex-valued.  It  has  been  assumed 
for  simplicity  that  x and  z are  of  the  same  dimension.  The 
final  result  is  obtained  by  applying  the  inverse  transforma- 
tion to  the  estimate  for  X computed  in  this  way.  It  can  be 
shown  (Ref.  26)  that  this  (unmodified)  estimate  agrees  pre- 
cisely with  the  one  that  would  have  resulted  had  no  transfor- 
mation been  interposed  at  all,  i.e.,  the  result  of  Eq . 2.2-3. 
The  objective  of  this  procedure  is  to  select  F to  reduce  the 
overall  cost  of  computing  estimates. 

The  numerical  structure  of  the  least-squares  equa- 
tions under  the  transformation  of  Eq.  3.2-1  determines  the 
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computational  cost  of  the  procedure.  The  transformed  least- 
squares  estimation  and  error  analysis  expressions  are  (Ref.  27) 
simply  the  usual  formulas,  Eqs.  2.2-3  and  2.2-4,  with  the 
transformed  statistics 


CXX  " FCxxF 


'ZZ 


FC  Ft 
zz 


'XZ 


FC  F' 
xz 


(3.2-2) 


An  appropriate  transformation  can  simplify  these  covariance 
matrices  and  the  resulting  equations  considerably.  It  is  well 
known  that  the  Karhunen-Loeve  transformation  (KLT)  is  optimal 
in  the  sense  that  it  diagonalizes  the  least-squares  equations 
(Ref.  28).  This  produces  the  greatest  reduction  in  the  cost 
of  computing  the  transformed  estimates.  Unfortunately,  this 
otherwise  ideal  transformation  is  very  expensive  to  compute. 
Recent  design  procedures  have  centered  on  the  use  of  alter- 
native transformations  which  combine  the  desirable  diagonali- 
zation  property  with  computational  efficiency. 

The  frequency  domain  estimation  procedures  of  this 
report  are  generalized  Wiener  filtering  algorithms  employing 
the  discrete  Fourier  transform.  This  transformation  can  be 
represented  by  the  n * n matrix  (Ref.  25,  p.  5) 

P . [Fjk]  =.  i[e-2!Tijk/nj  (3.2-3) 

where  0 <_  j,k  <.  n-1.  Two  properties  make  this  transformation 
attractive : 

• The  transformed  least-squares  equations 
resulting  from  Eq.  3.2-2  are  approxi- 
mately diagonal  for  a broad  class  of 
statistical  models. 


*The  superscript  (+)  denotes  the  complex  conjugate,  transpose 
operat ion. 
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• The  use  of  the  fast  Fourier  transform 
algorithm  allows  the  necessary  trans- 
formations to  be  carried  out  with  a high 
degree  of  efficiency. 

Frequency  domain  estimation  techniques  thus  draw  on  both 
statistical  theory  and  the  theory  of  "fast"  transforms  for 
their  utility.  The  statistical  theory  of  spectral  represen- 
tation, outlined  briefly  in  the  last  section,  provides  for  the 
approximate  diagonality  of  the  transformed  estimation  equations, 
while  the  FFT  technique  makes  the  transform  domain  efficiently 
accessible.  The  single-sensor  version  of  this  algorithm  is 
the  discrete  Wiener  filter  described  in  Section  3.1,  and 
presented  in  Fig.  3.2-1. 


OBSERVED  DATA 


H-24171 


ESTIMATES  AND 
STATISTICAL 
ERROR  ANALYSIS 


Figure  3.2-1  The  Generalized  Wiener  filter 


The  application  of  the  frequency  domain  estimation 
method  involves  neglecting  the  off-diagonal  terms  in  the  trans- 
formed covariance  matrices  of  Eq.  3.2-2.  These  terms  corre- 
spond to  cross  covariances  which  are  approximately  zero  accord- 
ing to  the  spectral  representation  theory  of  Section  3.1. 
Simplification  of  the  estimation  and  error  analysis  expressions 
by  neglecting  these  small  terms  is  responsible  for  the  sub- 
stantial savings  in  computer  time  associated  with  the  frequency 
domain  approach.  The  approximations  inherent  in  neglecting 
off-diagonal  terms  leads  to  what  have  been  called  "edge  effects" 
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(Ref.  32)  in  the  resulting  estimates.  A considerable  body  of 
matrix  theory  relates  to  the  analysis  of  these  approximations 
(Refs.  13,  14,  15  and  35).  The  analysis  briefly  outlined  in 
Appendix  A shows,  for  a given  estimation  problem  (of  the  sort 
which  obeys  the  restrictions  below),  that  edge  effects  become 
arbitrarily  small  for  a survey  track  of  sufficient  length. 

In  this  sense,  the  estimates  produced  using  frequency  domain 
techniques  are  asymptotically  equal  to  the  optimal,  least- 
squares  estimates.  The  impact  of  edge  effects  on  estimation 
error  is  investigated  numerically  in  Section  3.4. 

The  requirements  for  applicability  of  the  frequency 
domain  techniques  of  this  report  arise  from  three  sources: 

• The  need  for  regularity  conditions  to 
insure  strong  diagonality  in  the  fre- 
quency domain. 

• Appropriate  data  organization  for  the 
application  of  discrete  Fourier  trans- 
forms . 

• The  suitability  of  statistical  gravity 
models  which  are  stationary  in  a Cartesian 
coordinate  system. 

Note  that  this  last  condition  does  not  restrict  the  method 
to  "flat-earth"  geometry;  distance  along  the  curved  surface 
of  the  earth  can  be  used  as  the  independent  variable.  Very 
large  surveys  can  thus  be  accommodated.  A more  problem- 
oriented  list  of  restrictions  is  as  follows: 

• Statistical  models  for  the  relevant  grav- 
ity quantities  and  for  measurement  noise 
must  be  stat ionary . 

• Observations  must  be  taken  at  uniform 
spacing  (or  it  must  be  possible  to  mathe- 
matically place  the  problem  in  this  format). 
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• The  length  of  the  survey  track  must  be 

at  least  several  times  the  longest  corre- 
lation distance  appearing  in  the  statisti- 
cal models. 


3.3  NUMERICAL  SIMULATION 


An  estimation  error  covariance  analysis  is  described 
in  this  section  as  a numerical  demonstration  of  the  frequency 
domain  estimation  procedure.  The  problem  considered  is 
the  airborne  gradiometric  survey  of  gravity  anomalies  as 
described  in  Section  2.1.  The  rms  gravity  anomaly  estimation 
error  is  computed  as  a function  of  the  survey  parameters, 
including  gradiometer  accuracy  and  aircraft  altitude,  using 
Eq.  3.1-4.  This  numerical  work  is  carried  out  under  the 
following  assumptions: 


• For  simplicity,  only  a single  component 
of  the  gravity  gradient  tensor  is  ob- 
served. The  particular  component  used  is 
the  one  judged  to  be  most  closely  related 
to  gravity  anomaly. 

• The  required  covariance  function  models 
are  obtained  from  the  attenuated  white 
noise  model  described  in  Section  2.3. 

• Gradiometer  measurement  error  is  modeled 
as  white  noise. 

• The  observed  gravity  gradient  is  the  re- 
sult of  a ten  second  "moving  window" 
average  of  the  gradiometer  output. 


These  assumptions  are  discussed  below. 


The  measured  gravity  gradient  element  is  taken  to 
be  the  along-track  derivative  of  the  vertical  gravity  dis- 
turbance, defined  by 
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Y ( t ) = 5g(t,h)  (3.3-1) 

where 

6g  = vertical  component  of  the  gravity 
disturbance  vector 

t = distance  along  survey  track,  0 1 t 1 T 
h = aircraft  altitude 

This  seems  a natural  choice  because  of  the  close  relationship 
between  gravity  anomaly  and  vertical  gravity  disturbance,  and 
because  the  latter  quantity  would  be  determined,  within  a con- 
stant of  integration,  by  a precise  measurement  of  y(t)  along 
the  survey  track. 

The  required  gravity  model  for  this  problem  consists 
of  the  autocovariance  functions  for  the  anomalous  gravity 
gradient  at  the  aircraft  altitude  and  gravity  anomaly  at  the 
surface,  as  well  as  the  cross  covariance  function  between 
these  two  quantities.  These  covariance  functions  are  obtained 
from  the  attenuated  white  noise  model  using  automated  differ- 
entiation as  described  in  Section  2.3. 

The  modeling  of  gradiometer  measurement  error  as 
white  noise  is  consistent  with  the  error  models  that  have 
been  proposed  for  these  instruments  (Refs.  30,  31,  32).  This 
model  has  the  physical  appeal  that  it  is  a good  representation 
for  random  molecular  or  electronic  error,  which  is  recognized 
as  a limiting  error  source.  Gradiometer  measurement  accuracy 
is  specified  by  the  rms  error,  in  Eotvos  Units  (EU),  present 
in  the  moving  window  average  of  the  instrument  output.  The 
ten  second  averaging  time  was  chosen  to  be  in  agreement  with 
conventional  usage  in  the  gradiometer  community.  Reference  33 
contains  a thorough  treatment  of  the  relationships  among  white 
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noise  power  spectral  density,  gradiometer  averaging  time,  and 
reported  gradiometer  accuracy. 


Figure  3.3-1  displays  gravity  anomaly  estimation  error 
as  a function  of  gradiometer  accuracy  and  aircraft  altitude. 

The  numerical  value  of  twenty  thousand  feet  was  chosen  as  a 
representative  survey  altitude,  while  curves  for  surveys  on 
the  surface  and  at  forty  thousand  feet  are  shown  for  comparison. 
The  0-20  EU  accuracy  regime  displayed  includes  gradiometers  con- 
sidered to  be  of  moderate  accuracy,  as  well  as  ideal,  noise- 
free  gradiometry.  Note  that  the  data  for  this  simulation 
consists  of  gradiometry  only;  no  initial  gravity  anomaly 
"fix"  has  been  assumed.  It  is  recognized  that  practical 
gravity  surveys  would  involve  some  direct  measurement  of  the 
gravity  anomaly.  An  interesting  feature  of  these  results  is 
that  they  show  a significant  residual  estimation  error  with 
a noise-free  gradiometer.  This  can  be  explained  by  the  fact 
that  the  gradiometer,  as  a differentiating  device,  cannot 
measure  the  bias  (constant)  portion  of  the  gravity  anomaly. 

For  finite  track  lengths  this  inability  to  measure  the  zero 
frequency  component  also  extends  to  the  very  low  frequencies. 

The  estimation  error  power  spectral  density  function 
describes  the  "frequency  response"  of  the  gradiometer  as  a 
device  for  estimating  gravity  anomaly.  As  discussed  in 
C^ction  3.1,  the  discrete  power  spectral  density  4>ee(w), 
evaluated  at  wave  number  u),  can  be  viewed  as  the  variance  of 
the  component  of  the  estimation  error  at  a given  frequency. 

The  theory  of  spectral  representation  also  states  that  the 
components  of  the  estimation  error  at  any  two  different 
frequencies  are  statistically  uncorrelated.  This  is  consis- 
tent with  Parseval's  identity,  Eq.  3.1-5,  which  is  simply  the 
statement  that  the  total  mean  square  estimation  error  is  the 
sum,  across  all  wave  numbers  - of  $ee(u)'  These  characteristics 
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Figure  3.3-1  Projected  Gravity  Anomaly  Estimation  Error 

for  an  Airborne  Survey  with  Gradiometry  Only 


make  ♦ (<*>)  a useful  representation  of  the  estimal  ion 

error . 


Figure  3.3-2  compares  the  estimation  error  power 
spectral  density  function  for  1 and  5 EU  gradiometer 
accuracies  with  the  gravity  anomaly  power  spectral  density 
from  which  it  is  derived.  As  anticipated,  the  estimation 
error  spectrum  is  almost  entirely  confined  to  the  very  low 
wave  numbers,  reflecting  the  gradiometer ' s insensitivity  to 
the  bias  component  of  the  gravity  anomaly.  This  behavior 
is  seen  even  more  clearly  in  the  spectral  coherency  func- 
tion (Ref.  21)  defined  in  the  notation  of  Section  3.1,  by 


2 

P 


(<D) 


$ ( 0)  ) 
*xz 


*xx(u,)  *zz(ul) 


(3.3-2) 
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WAVE  NUMBER 

(NUMBER  OF  WAVELENGTHS  PER  TRACK  LENGTH! 


Figure  3.3-2  Power  Spectra  for  Airborne  Gradiometry 

■ 

I 


The  usefulness  of  the  spectral  coherency  function  is  due  to 
the  fact  that  Eq . 3.1-4  may  be  rewritten  as 


*ee(u,)  = f1  ‘ p2(w>]  *xx<“>  (3.3-3) 


a 

i 


i 


j 


Thus  the  spectral  coherency  function  defines  a kind  of  trans- 
fer function  between  the  gravity  anomaly  power  spectrum  and 
the  estimation  error  spectrum  of  the  gravity  anomaly.  A plot 
of  the  gravi ty-anomaly/gravity-gradient  spectral  coherency 
function  for  rms  gradiometer  accuracies  of  1 and  5 EU  is 
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presented  in  Fig.  3.3-3.  These  curves  show  the  effect  of 
increased  gradiometer  noise  as  increased  "roll-off"  at  very 
low  and  moderately  high  wave  numbers.  The  impact  this  roll 
off  has  on  the  estimation  error  spectrum  is  weighted  by  the 
gravity  anomaly  power  spectrum  (Eq.  3.3-3).  The  high  freq- 
uency roll-off  has  little  effect,  as  is  apparent  from  Fig. 
3.3-2,  because  the  gravity  anomaly  has  very  little  energy 
in  this  part  of  the  spectrum.  However,  the  low  frequency 
roll-off  in  spectral  coherency  has  a strong  effect  on  esti- 
mation error  because  the  gravity  anomaly  has  a substantial 
low  frequency  component.  The  spectral  coherency  function 
provides  a means  of  assigning  changes  in  the  estimation 
error  spectrum  directly  to  changes  in  sensor  response,  as 
weighted  by  the  power  spectrum  of  the  estimated  variable. 

RMS  GRADIOMETER  ERROR 

1.0  | ~ :■  ~ — — l EU 


\ 6 EU 


• RMS  GRAVITY  ANOMALY  - Si 

• ANOMALY  CORRELATION  DISTANCE  - 30  nm 

• GRAOIOMETER  AVERAGING  TIME  • 10  >*c 

• AIRCRAFT  SPEED  • SOD  ki 

• AIRCRAFT  ALTITUDE  - 30,000  It 

• SURVEY  TRACK  LENGTH  • 1300  nm 


WAVE  NUMBER 

(NUMSER  OP  WAVELENGTHS  PER  TRACK  LENGTH! 


Figure  3.3-3 


Gravi ty -Anomaly /Gravi ty-^radient 
Spectral  Coherency  Functions  for 
Airborne  Gradiometry 
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3.4  EDGE  EFFECT  ANALYSIS 

For  a given  estimation  problem  the  gravity  quantity 
estimates  commputed  with  frequency  domain  methods  may  differ 
from  those  obtained  via  direct  least-squares  estimation.  From 
the  spectral  theory  viewpoint  presented  in  Section  3.1,  this 
difference  is  because  the  Wiener-Kolmogorov  filter  is  optimal 
only  for  the  case  of  a survey  track  of  infinite  length.  Edge 
effect  approximations  result  when  the  filter  is  used  for  sur- 
vey tracks  of  finite  length.  In  the  generalized  Wiener  filter 
view,  edge  effects  are  the  result  of  neglecting  off-diagonal 
terms  in  the  transformed  estimation  equations.  From  this 
perspective,  edge  effects  are  the  result  of  a tradeoff  of 
increased  estimation  error  for  reduced  processing  time.  Of 
course,  the  objective  in  designing  such  an  estimation  proced- 
ure is  to  obtain  a very  substantial  reduction  in  processing 
time  in  exchange  for  a very  small  increase  in  estimation  error. 

The  analysis  outlined  in  Appendix  A shows  that  the 
increase  in  rms  estimation  error  due  to  edge  effects  becomes 
arbitrarily  small  as  survey  track  length  is  increased.  This 
result  provides  theoretical  justification  for  the  use  of  freq- 
uency domain  techniques.  The  purpose  of  this  section  is  to 
present  numerical  verification  of  this  result  for  a survey 
track  of  finite  length. 

Numerical  evaluation  of  the  portion  of  estimation 
error  due  to  edge  effects  can  be  obtained  by  exercising  the 
appropriate  error  analysis  expression.  The  estimation  error 
covariance  matrix  associated  with  an  arbitrary  linear  esti- 
mation procedure  of  the  form 

x = Gz  ( 3 . 4-1 ) 
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is  given  by  (Ref.  17) 


(3.4-2) 


The  least-squares  error  covariance  expression,  Eq.  2.2-4, 
is  a special  case  of  Eq . 3.4-2,  when  G is  the  optimal  gain 
matrix  of  Eq.  2.2-3.  With  the  appropriate  choice  for  G, 

Eq.  3.4-2  describes  the  expected  squared  error  associated 
with  frequency  domain  estimation  including  edge  effect 
approximations.  The  proper  G matrix  is  obtained  through  the 
procedure  used  to  construct  the  frequency  domain  filter  as 
block-diagrammed  in  Fig.  3.2-1. 


In  Fig.  3.4-1  the  rms  gravity  anomaly  estimation 
error  including  edge  effects  (Eq.  3.4-2)  is  compared  with  the 
optimal  estimation  error  curve  from  Fig.  3.3-1.  Note  that 
the  curve  labeled  "optimal  estimation"  does  not  include  the 
edge  effect  degradation.  In  the  worst  case,  the  increase  in 
estimation  error  over  the  optimum  is  a small  fraction  of  the 
modeled  rms  gravity  anomaly  level . The  increased  estimation 
error  resulting  from  edge  effects  becomes  more  prominent  as 
gradiometer  accuracy  is  increased.  This  may  be  interpreted 
as  a result  of  the  greater  demands  for  precision  placed  on 
the  estimation  procedure  by  the  availability  of  highly  accu- 
rate data. 
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4.  MULTISENSOR  FREQUENCY  DOMAIN  ESTIMATION 

4.1  PROBLEM  FORMULATION 

Frequency  domain  estimation  for  multisensor  problems 
involves  some  data  structure  considerations  not  present  in 
the  single  sensor  case.  Mathematical  structure  in  the  least- 
squares  estimation  equations  appropriate  for  the  use  of  fre- 
quency domain  methods  depends  upon  relationships  between  the 
data  grids  for  the  various  sensors.  The  treatment  of  multi- 
sensor frequency  domain  estimation  is  begun  in  this  section 
by  dividing  these  problems  into  two  classes  according  to  the 
way  the  data  has  been  taken.  For  one  class  of  problems, 
where  the  sensor  data  grids  are  all  "matched"  in  a certain 
sense,  the  applicable  frequency  domain  algorithm  does  not 
differ  in  any  fundamental  way  from  the  single-sensor  method. 
The  necessary  extension  of  the  techniques  of  Chapter  3 is 
straightforward,  with  detailed  formulas  presented  in 
Appendix  B. 

The  second  class  of  problems,  characterized  by  sensor 
data  grids  which  are  "mismatched",  is  more  difficult,  and  its 
solution  becomes  the  major  theme  of  the  next  two  sections. 
Section  4.2  presents  a sequential,  sensor-by-sensor  estima- 
tion procedure.  While  this  technique  by  itself  does  not 
result  in  a faster  algorithm,  it  does  provide  a useful  frame- 
work around  which  to  organize  multisensor  procedures.  The 
derivation  of  frequency  domain  statistics  suitable  for  use  in 
multisensor  estimation  is  outlined  in  Section  4.3.  The  use 
of  a discrete  version  of  the  classical  Fourier  transform 
sampling  theorem  is  the  major  part  of  this  development.  The 
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application  of  this  result  in  conjunction  with  the  sequential 
estimation  approach  provides  a solution  to  the  second,  more 
difficult  class  of  multisensor  estimation  problems. 

Throughout  this  report  it  is  assumed  that  the  data 
are  observed  on  one-dimensional  grids  with  uniform  spacing, 
or,  at  least,  that  the  situation  appears  this  way  in  terms  of 
the  mathematics  of  estimation.  In  treating  multisensor 
estimation,  it  is  convenient  to  distinguish  between  problems 
where  the  data  grids  for  the  various  sensors  are  matched, 
and  where  they  are  not.  In  this  context,  two  data  grids  are 
matched  if  they  have  these  characteristics: 

• Both  grids  have  uniform  measurement 
point  spacing.  The  point  spacing  need 
not  agree  between  grids. 

• Both  grids  have  the  same  number  of  data 
points. 

This  definition  clearly  includes  grids  which  are  physically 
matched,  i.e.,  those  for  which  the  measurement  points  fron 
different  grids  overlay  one  another.  It  also  encompasses 
those  situations  where  data  grids  are  matched  mathemat ical ly , 
perhaps  involving  a change  of  scale  or  a phase  shift. 

Gravity  gradiometry  involving  the  measurement  of 
more  than  one  gradient  provides  a ready  example  of  the  matched 
grid  estimation  problem.  When  a moving-base  gradiometer  is 
used  to  make  a series  of  simultaneous  measurements  of  two  or 
more  gradient  elements,  the  data  from  each  of  these  sensors 
fall  at  the  same  spatial  points.  The  resulting  data  grids 
are  physically  matched. 
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The  use  of  airborne  gradiometry  in  conjunction  with 
point  gravimetry  at  the  surface  illustrates  the  estimation 
problem  with  mismatched  data  grids.  Airborne  gradiometer 
measurements  are  available  on  a continuous,  or  nearly  con- 
tinuous, basis.  However,  surface  gravimetry  points  may  be 
separated  by  many  miles,  or  many  tens  of  miles.  For  a sur- 
vey of  the  same  physical  region,  the  data  grids  from  these 
two  sensors  cannot  be  made  compatible.  Multisensor  surveys 
with  mismatched  data  grids,  such  as  this  one,  constitute 
an  important  class  of  estimation  problems.  Frequency 
domain  techniques  for  dealing  with  these  problems  are 
developed  in  the  next  two  sections. 

The  single-sensor  frequency  domain  methods  of 
Chapter  3 are  readily  adapted  to  multisensor  problems  of  the 
sort  with  matched  data  grids.  This  can  be  accomplished  by 
treating  the  observed  data  as  a vector-valued,  rather  than  a 
scalar,  function  of  distance  along  the  survey  track.  The 
generalization  in  the  frequency  domain  estimation  expres- 
sions then  simply  involves  the  replacement  of  certain  scalar 
quantities  with  their  vector  or  matrix  analogs.  This  exten- 
sion is  made  explicit  in  Appendix  B.  An  equivalent  estima- 
tion procedure,  based  on  sequential  techniques,  is  discussed 
in  the  next  section. 

4.2  SEQUENTIAL,  SENSOR-BY-SENSOR  ESTIMATION 

Sequential  data  processing,  handling  data  from  only 
one  sensor  at  a time,  helps  to  mitigate  the  added  complexity 
of  multisensor  estimation.  The  sequential  estimation  and 
error  analysis  procedure  presented  in  this  section  is  well 
known  and  has  found  application  in  many  forms  (Refs.  17,  29, 
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35).  It  should  be  stressed  that  the  technique  explicitly  con- 
sidered here,  unlike  the  recursive  Kalman  filter  algorithm, 
does  not  directly  reduce  computational  effort.  Its  main 
utility  is  rather  that  it  reduces  the  complexity  of  the  esti- 
mation and  error  analysis  expressions  that  must  be  considered 
at  any  one  time.  In  conjunction  with  frequency  domain  method- 
ology, this  is  advantageous  because  it  helps  to  expose  the 
mathematical  structure  which  accompanies  each  sensor. 

Consider  the  estimation  of  an  unknown  vector  x using 
the  information  contained  in  two  data  vectors  u and  v.  The 
unknown  vector  might  consist  of  gravity  anomaly  values  on  some 
survey  grid,  while  u and  v contain  the  data  from  two  different 
sensors  such  as  gravity  gradiometry  and  gravimetry.  The  least- 
squares  estimate  for  x based  only  on  the  information  in  u and 
the  accompanying  estimation  error  covariance  matrix  are,  from 
Eqs.  2.2-3  and  2.2-4: 


x = C C-1  u 
— XU  uu  — 


,-l 


Cee  Cxx  ” Cxu  Cuu  Cux 


(4.2-1) 

(4.2-2) 


The  results  of  Eqs.  4.2-1  and  4.2-2  can  be  updated  recursively 
to  account  for  the  additional  information  provided  by  v.  If 
xu  and  Cee<u  are  used  to  denote  the  results  of  Eqs.  4.2-1  and 
4.2-2,  respectively,  then  the  statistical  literature  (Ref.  35) 
provides  the  estimate  and  error  matrix  based  on  both  u and  v 
as 


— uv  — u 


,-l 


,-l 


= x..  + C C (v  - c C u) 

XV  *U  wuv—  VU  Ulh ' 


,-l 


c = c -c  c C 

ee*uv  ee*u  xv*u  vv*u  vx*u 


(4.2-3) 

(4.2-4) 


4-4 


THE  ANALYTIC  SCIENCES  CORPORATION 


These  equations  represent  another  application  of  Eqs.  2.2-3 
and  2.2-4,  except  that  the  covariance  matrices  involved  are 
conditioned  on  u.  That  is,  each  covariance  matrix  is  present 
in  its  usual  context,  with  its  form  modified  to  account  for 
the  information  in  u.  These  conditional  covariance  matrices 
are  given  by  the  expressions 


C = C 

vvu  vv 

C = C 

XV'U  xv 


C C"1  C 
vu  uu  uv 


C C"1  C 
XU  uu  uv 


(4.2-5) 

(4.2-6) 


Note  that  the  single  sensor  error  covariance  expression  (Eq. 

4.2- 2)  is  also  the  conditional  covariance  expression  for  C 

ee  • u 

as  might  be  expected.  Reference  36  discusses  the  statistical 
interpretation  of  covariances  conditioned  on  the  data  already 
processed . 

The  procedure  described  above  extends  naturally  from 
two  to  an  arbitrary  number  of  sensors.  The  use  of  conditional 
covariance  expressions  analogous  to  Eqs.  4.2-2,  4.2-5  and 

4.2- 6  provides  conditional  statistics  based  on  the  group  of 
sensors  already  processed.  The  substitution  of  the  resulting 
covariance  matrices  into  the  conventional  estimation  formulas 
Eqs.  2.2-3  and  2.2-4  yields  the  estimation  and  error  analysis 
expressions  for  incorporating  the  next  sensor.  If  one  chooses 
to  view  the  vector  u as  containing  the  data  from  all  the  sen- 
sors processed  up  to  a given  point,  with  v constituting  the 
data  from  the  next  sensor  to  be  added,  then  Eqs.  4.2-3  and 

4.2- 4  are  the  general  iterative  expressions  just  described. 
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4.3  SENSOR-BY-SENSOR  ESTIMATION  IN  THE  FREQUENCY  DOMAIN 

Consider  the  two  sensor  estimation  problem  of  the 
last  section,  where  x is  to  be  estimated  using  data  contained 
in  vectors  u and  v.  First  assume  for  simplicity  that  x,  u, 
and  v are  all  of  dimension  n.  Then  the  finite  Fourier  trans- 
form of  Eq.  3.2-3  may  be  applied  to  each  of  these  vectors, 
producing  the  transforms 

X = Fx,  U = Fu,  V = Fv  (4.3-1) 

According  to  the  generalized  Wiener  filter  theory  of  Section 
3.2,  the  covariance  matrices  associated  with  these  transformed 
quantities  may  be  taken  to  be  diagonal,  with  diagonal  terms 
corresponding  to  the  appropriate  discrete  power  spectral 
density.  From  Section  3.2,  these  matrices  are  all  of  the 
form* 


CXX  = F CxxF+  = diag(*xx(0) *xx(n-l))  (4.3-2) 


The  substitution  of  the  diagonal  covariance  matrices 
into  the  sensor-by-sensor  estimation  equations,  Eqs.  4.2-1 
through  4.2-4,  yield  the  set  of  equivalent  scalar  frequency 
domain  expressions  below. 


X (u>)  = $xu(w)  U(u>) 

<t>  (w) 

uuv  ’ 


ee*u 


(w) 


*xx<"> 


l*xu(l*» 

*uu(w> 


(4.3-3) 


(4.3-4) 


♦The  notation  A = diar( ai , . . . , an ) indicates  the  matrix  with 
elements  Aii  = at , Aij  = 0 ; 1 <_  i,j  £ n,  i^j . 


I 


J 
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X (u)  = X <«>)  + *xvu(M)  [V(u)  - 4'vu(l'l)  <4.3-5) 

UV  U » TT37 


vvu 


uu 


4>  (u)  = $ (u>)  - 

ee*uvv  ee*u 


3>  (u)r 

XV«Uv  71 

<f>  (ui) 

vvu  ’ 


(4.3-6) 


A similar  substitution  into  the  conditional  covariance 
expressions,  Eqs.  4.2-5  and  4.2-6  provides 


4>  (oj)  = 
vvu 


$ ( Ul  ) 

vv'  ' 


*uu<"> 


$ (w) 

xvu 


<J>  ( U)) 

xv'  ' 


$ (u)4>  (u)) 

XUV  UV 

♦uu'"* 


(4.3-7) 


The  use  of  Eqs.  4.3-3  through  4.3-7  in  an  integrated, 
muitisensor  algorithm  is  presented  in  Fig.  4.3-1.  The  proce- 
dure of  the  figure  provides  estimates  and  estimation  error 
analysis  based  on  an  arbitrary  number  of  sensors,  with  the 
restriction  that  all  of  the  sensor  grids  must  be  of  the 
matched  type.  This  recursive  algorithm  is  mathematically 
equivalent  to  the  batch  procedure  of  Appendix  C. 

The  recursive  algorithm  of  Fig.  4.3-1  can  be  extended 
to  many  types  of  problems  with  mismatched  sensor  data  grids  by 
applying  a sampling  theorem.  This  is  used  to  deduce  the  for- 
mat of  the  appropriate  frequency  domain  expressions  for  use  in 
the  sequential  procedure  of  Section  4.2.  Let  v and  u be  vec- 
tors of  dimension  m-n,  representing  sensor  data  sampled  on 
mismatched  grids  such  that  n = rm,  where  r is  an  integer. 

The  implications  of  this  assumption  for  practical  estimation 
are  discussed  below.  Under  these  circumstances  v may  be 


4-7 


the  analytic  sciences  corporation 


regarded  as  a sampling  taken  from  an  n-vector  v*  with  an 
artificial  "data  grid"  matched  to  that  of  u.  The  elements 
of  v are  related  to  those  of  v*  by 


v(j) 


v* 


( j- 


1 )r  + 


j=l,  . . . ,m 


(4.3-8) 


U 


<W  Cxu'Cxx 


M42U 


V 


c c c 

UV  XV'  W 


XUM 


Figure  4.3-1  Rensor-by-Sensor  Frequency 

Domain  Estimation 

This  approach  has  the  advantage,  according  to  the  spectral 
theory  of  Chapter  3,  that  the  frequency  domain  covariance 
matrices  involving  v*  are  of  known  diagonal  form.  The  corre- 
sponding matrices  for  v can  then  be  computed  directly  from 
the  algebraic  relationship  between  v and  v*  in  the  frequency 
domain . 
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(Ref.  25)  provides  the  necessary  connection  between  the  fre- 
quency domain  versions  of  v and  v*.  Letting  V(w),  io=0 

m-1  and  V*(w),  uj=0 n-1  denote  these  transforms,  the  samp- 

ling theorem  provides 

r-1 

V(u>)  * l V*(oj+mp)  (4.3-9) 

p=0 

Figure  4.3-2  illustrates  the  sampling  theorem  relationship. 
The  algebra  of  the  application  of  Eq.  4.3-9  to  the  calcula- 
tion of  the  statistics  of  V(w)  is  treated  in  Appendix  C. 


V*  (j)  v'<«) 


Figure  4.3-2  Sampling  Theorem  for  Finite  Transforms 


The  structure  of  the  multisensor  estimation  algori- 
thm presented  in  Fig.  4.3-1  remains  appropriate  for  problems 
with  mismatched  data  grids.  Certain  of  the  spectral  equations 
referred  to  in  the  "SENSOR  ADDITION"  block  must  be  replaced 
to  account  for  the  sampled  data  analysis  discussed  above.  From 
Appendix  C,  the  spectral  estimation  and  error  analysis  formulas 
analogous  to  Eqs.  4.3-5  and  4.3-6  are 
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xuv(u)  ' Vw>  + *xv*(u)'*vv-u(“  mod  n)‘(V-%u*uiU)(u)  mod  n) 


2 —1 

4>  (iii)  = <t » (id)-[<J>  ^ ( <d  ) I • 4>  (uj  mod  n ) 

ee*uv  ’ ee*u  1 xv**uv'  ’ 1 vv*u  ' 


(4.3-10) 


(4.3-11) 


for  u)=0 n-1.  The  expression  (w  mod  m)  denotes  the  remain- 

der after  integer  division  by  m, 

i.e.,  (jm+k)  mod  m = k 

for  any  integer  j and  any  0 < k < m-1. 


The  conditional  spectral  densities  which  replace  Eq . 
4.3-7  for  the  mismatched  sensor  grid  case  are  given  by 


*wu(to)  = 2 $v*v*(03+mp)  - 

p=0 


*uv*(“+mp) I * 

$>uu(u)+mp) 


(4.3-12) 


for  0 £ oj  <.  m-1,  and 


<txv*-u(a))  = 4>xv*(ui)  ‘ <txu(li,)$uu(ll>Huv*(u)) 


(4.3-13) 


for  0 <.  id  <.  n-1.  The  substitution  of  Eqs.  4.3-10  and  4.3-11 
for  Eqs.  4.3-5  and  4.3-6;  and  the  replacement  of  Eq.  4.3-7  by 
Eqs.  4.3-12  and  4.3-13  in  the  procedure  of  Fig.  4.3-1  allow 
for  the  processing  of  mismatched  multisensor  data. 

The  multisensor  estimation  algorithms  described  in 
this  section  are  subject  to  two  restrictions.  One  of  these 
has  been  applied  to  simplify  the  exposition,  the  other  is  more 
essential.  It  has  been  assumed  that  the  vectors  x and  u are 
of  the  same  dimension.  This  is  equivalent  to  assuming  that 
these  quantities  have  been  sampled  on  matched  grids,  which 
serves  to  simplify  the  exposition.  This  restriction  may 
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be  removed  by  applying  the  sampling  theorem  approach  for 
mismatched  grids  which  has  just  been  described.  A second, 
more  basic  restriction  on  the  algorithm  is  the  assumption, 
in  the  mismatched  data  grid  case,  that  the  ratio  of  the 
amounts  of  data  corresponding  to  any  two  sensor  grids  (i.e. 
the  number  r in  the  analysis  above)  is  an  integer.  In  many 
problems  of  interest  this  restriction  can  be  met  through  an 
appropriate  arrangement  of  the  input  data.  For  example,  when 
a data  "mismatch"  results  from  discretizing  an  essentially 
continuous  measurement,  such  as  airborne  gradiometry,  pre- 
processing may  be  applied  to  alter  the  effective  sampling 
rate.  This  may  be  accomplished  with  no  information  loss  by 
smoothing  with  a moving  window  averager,  and  then  sampling 
at  a carefully  chosen  increased  rate.  The  new  sampling  rate 
may  be  chosen  so  that  the  data  mismatch  ratio,  r,  is  an 
integer.  Frequency  domain  techniques  are  then  directly 
applicable.  With  suitable  structuring  of  the  input  data,  the 
computational  efficiency  associated  with  frequency  domain 
methodology  is  available  over  a wide  range  of  multisensor 
estimation  problems. 


4.4  MULTI SENSOR  NUMERICAL  SIMULATION 

The  results  of  numerical  simulations  with  the  multi- 
sensor frequency  domain  estimation  procedure  are  presented 
in  this  section  to  demonstrate  the  applicability  of  the  method. 
These  results  consist  of  frequency  domain  estimation  error 
analyses  of  two-  and  three-sensor  simulated  gravity  anomaly 
surveys.  Simulated  data  sources  for  the  two-sensor  survey 
consist  of  airborne  gradiometry  and  surface  gravimetry.  The 
three-sensor  survey  is  formed  by  the  addition  of  satellite 
altimetry.  The  relevant  survey  geometry  for  these  simulations 
is  depicted  in  Fig.  1.2-1. 
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Apart  from  the  use  of  multisensor  techniques,  the 
methodology  applied  in  this  section  is  much  the  same  as  that 
described  in  Section  3.3.  The  statistical  models  and  assump- 
tions associated  with  the  simulation  of  airborne  gradiometry 
are  identical  with  those  of  that  section.  Similar  use  is 
made  of  the  attenuated  white  noise  model  in  the  automated 
generation  of  the  self-consistent  disturbance  gravity  models 
required  in  the  simulation.  This  procedure  has  also  been 
described  in  Section  2.3. 

Surface  gravimetry  is  modeled  as  a series  of  point 
gravity  anomaly  measurements  parameterized  by  the  (uniform) 
spacing  between  observation  points.  A nominal  rms  gravimetry 
accuracy  of  1 mgal  is  assumed  throughout.  Gravimetry  errors 
are  assumed  to  be  uncorrelated  between  observation  points. 

Satellite  radar  altimetry  observations  are  simulated 
as  noisy  measurements  of  the  undulation  of  the  geoid  beneath 
the  satellite's  orbit.  The  separation  between  successive 
undulation  measurements  is  parameterized  by  the  altimeter 
pulse  repetition  rate  and  the  satellite  orbital  velocity.  An 
orbital  velocity  of  8 km/sec  is  assumed  throughout.  Altimeter 
accuracy,  as  reflected  in  the  rms  accuracy  of  the  geoid  height 
determination,  is  a parameter  in  the  simulation. 

The  results  of  the  gradiometric/gravimetric  survey 
simulation  are  presented  in  Fig.  4.4-1.  An  interesting  feature 
of  this  figure  is  the  reduction  in  gravity  anomaly  estimation 
error  resulting  from  the  introduction  of  relatively  few  gravi- 
metry observations.  The  spectral  analysis  of  gradiometer  re- 
sponse presented  in  Section  3.3  has  shown  that  most  of  the 
estimation  error  after  a survey  with  a gradiometer  alone  is 
at  the  bias,  or  very  low  frequency,  level.  Gravimetry  data 
complements  gradiometry  very  well  by  providing  measurements 
sensitive  to  the  bias  component  in  the  gravity  anomaly. 
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KM*  GRAVITY  ANOMALY  • B4  mfrt 

ANOMALY  CORRELATION  DISTANCE  • 20  nm 

FOCT-REOUCTION  GRAVIMETRY  ERROR  • t meal  (WHITE! 

GRADIOMETER  AVERAGING  TIME  • 10  Me 

AIRCRAFT  SPEED  - 500  kt 

AIRCRAFT  ALTITUDE  • 20,000  It 

SURVEY  TRACK  LENGTH  • 1200  nm 


Figure  4.4-1  Projected  Gravity  Anomaly  Estimation  Error  for 

a Combined  Gradiometric/Gravimetric  Survey 

The  projected  gravity  anomaly  estimation  error  re- 
sulting from  two-  and  three-sensor  surveys  involving  gravi- 
metry, gradiometry,  and  satellite  altimetry  in  various  com- 
binations are  shown  in  Table  4.4-1.  The  two  sections  of  this 
table  contrast  the  quality  of  the  estimates  that  can  be  ob- 
tained using  a high  quality  gradiometer  (1  EU)  with  those 
obtainable  using  a gradiometer  of  lesser  accuracy  (10  EU). 

It  is  recognized  that  the  rms  estimation  errors  pre- 
sented in  Table  4.4-1  are  well  in  excess  of  current  gravity 
mapping  needs  (Ref.  38).  These  results  indicate  that  gradiom- 
etery  data  from  a single,  one-dimensional  track  is  not  adequate 
for  high  precision  gravity  anomaly  estimation.  The  addition 
of  multiple,  parallel  data  tracks  over  the  survey  area  should 
produce  considerable  improvement  in  gradiometer-based  gravity 
estimates. 
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TABLE  4.4-1 

PROJECTED  MULTISENSOR  GRAVITY  SURVEY  ACCURACIES 


SURVEY  WITH  HIGHLY  ACCURATE  GRADIOMETER 


■ ONE-DIMENSIONAL  SURVEY  GEOMETRY 

• RMS  GRAVITY  ANOMALY  - 64  mg«l 

■ ANOMALY  CORRELATION  DISTANCE  - 20  nm 

■ POST-REDUCTION  GRAVIMETRY  ERROR  - 1 mgal  (WHITE) 

• GRADIOMETER  AVERAGING  TIME  - 10«#c 

• AIRCRAFT  SPEED  - 600  kt 

• AIRCRAFT  ALTITUDE  - 20,000  ft 

■ SATELLITE  ALTIMETER  DATA  RATE  - 10/we 

• SURVEY  TRACK  LENGTH  • 1200  nm 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS 

* 

This  report  describes  the  development  of  new, 

f*  i • 

computationally-ef f icient  techniques  for  processing  multi- 
sensor gravity  data.  These  procedures  apply  modern  frequency 
domain  signal  processing  methodology  to  existing  statistical 
least-squares  gravity  estimation  techniques.  The  general 

i;  ' 

statistical  least-squares  approach  to  gravity  quantity 
estimation : 

• Produces  theoretically  optimum  gravity 

quantity  estimates 

• Applies  to  a very  general  class  of 
gravity  sensors  and  survey  types 

• Contains  a "built-in"  capability  for 

statistical  estimation  error  analysis 

• Requires  self-consistent  statistical 

gravity  models. 

The  conventional  least-squares  procedure  has  the  disadvantage 
that  numerical  processing  requirements  place  important 
limits  on  its  practical  utility.  A direct  implementation 
of  the  least-squares  equations  leads  to  computer  processing 
requirements,  as  well  as  susceptability  to  numerical  in- 
accuracy, that  sharply  restrict  problem  size.  However,  the 
new  procedure  largely  eliminates  these  practical  difficulties 
through  the  introduction  of  computationally  efficient  methods. 

The  frequency  domain  estimation  techniques  described 
in  this  report  are  designed  to  substantially  reduce  the  com- 
puter processing  requirements  of  the  statistical  least-squares 
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procedure  for  a wide  variety  of  gravity  applications.  This 
technique  is  based  upon  the  transformation  of  the  least-squares 
equations  into  a form  advantageous  for  computation.  The 
overall  cost  of  the  resulting  procedure  is  proportional  to 
the  cost  of  a fast  Fourier  transform  (FFT)  of  the  input  data. 

As  a result,  the  computer  burden  imposed  by  the  frequency  do- 
main method  is  less  by  orders  of  magnitude  than  that  of  the 
equivalent  direct  least-squares  procedure. 

The  efficiency  of  the  frequency  domain  method  is 
based  on  the  exploitation  of  the  mathematical  structure 
ordinarily  present  in  gravity  quantity  estimation  problems. 
Requirements  for  the  applicability  of  the  method,  based  on 
preserving  appropriate  structure  in  the  estimation  equations, 
are: 

• Regularly  gridded  data.  Data  may  be 
preprocessed  by  averaging  or  inter- 
polation. 

• Stationary  statistical  gravity  models. 

• Long  gravity  survey  track  length  compared 

to  the  relevant  gravity  quantity  correlation 
distances . 

The  class  of  problems  considered  in  this  report  has  been 
limited  to  estimation  along  and  above  a one-dimensional  sur- 
vey track.  It  is  anticipated  that  the  frequency  domain  tech- 
nique can  be  extended  for  use  with  two-dimensional  survey 
regions  under  requirements  similar  to  those  above. 

Numerical  simulations  have  been  carried  out  to 
demonstrate  the  performance  of  the  method.  Numerical  results 
include : 
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• Estimation  error  analyses  for  gravity 
anomaly  surveys  with  one,  two,  and  three 
sensors . 

• Preliminary  analyses  of  the  spectral  re- 
sponse characteristic  of  airborne  gravity 
gradiometry . 

• Investigation  of  estimation  error  sensi- 
tivity to  frequency  domain  edge  effect 
approximat ions . 

• Verification  of  estimated  computer  process- 
ing times  for  the  frequency  domain  method. 


The  results  presented  in  this  report  represent  a 
first  step  in  the  application  of  modern  frequency  domain  sig- 
nal processing  methodology  to  the  processing  of  gravity  data. 

A number  of  areas  suggest  themselves  as  being  worthy  of  future 
investigation : 


• Development  of  frequency  domain  methods 
suitable  for  the  estimation  of  gravity  on 
a two-dimensional  surface. 

• Investigation  of  frequency  domain  data  com- 
pression techniques.  The  results  obtained 
with  frequency  domain  least  squares  methods 
suggest  a similar  approach  to  the  efficient 
compression  of  gravity  data. 

• Investigation  of  edge-effect  compensation 
procedures.  Numerical  results  on  estima- 
tion error  sensitivity  to  the  effects  of 
finite  track  length  (Section  3.4)  suggests 
the  use  of  modified  frequency  domain  tech- 
niques for  some  problems. 

• Development  of  spectral  sensor  response 
characterization  techniques  appropriate 
for  multisensor  applications.  Sensor 
response  comparisons,  analysis  of  sensor 
information  overlap,  and  multisensor  trade- 
off studies  should  benefit  from  an  exten- 
sion of  the  coherency  function  idea  applied 
to  gradiometry  in  Section  3.4. 
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APPENDIX  A 

ASYMPTOTIC  OPTIMALITY  OF  FINITE, 


FREQUENCY  DOMAIN  ESTIMATION 


This  appendix  presents  a brief  analysis  of  the  asymp- 
totic optimality  property  of  the  finite,  frequency  domain 
estimation  methods  developed  in  this  report.  Under  appro- 
priate regularity  conditions,  it  is  shown  that  the  rms  est ima- 
tion  error  associated  with  the  frequency  domain  algorithm 
converges  to  the  optimum  as  the  survey  track  length  is  in- 
creased. This  analysis,  which  does  not  appear  to  be  available 
in  the  open  literature,  provides  explicit,  theoretical  justi- 
fication for  the  use  of  frequency  domain  procedures  in  con- 
nection with  finite  estimation  problems.  The  estimation  of 
gravimetric  quantities  from  data  is  such  a problem. 

For  the  problems  considered  in  this  report,  the  co- 
variance  matrices  appearing  in  the  estimation  equations  belong 
to  a class  of  matrices  known  as  Toeplitz  forms  (Ref.  37).  An 
nxn  matrix  C is  a Toeplitz  matrix  if  there  is  some  function 


♦ CD, 


-(n-1)  < Z < (n-1 ) 


( A-l ) 


for  which  the  (j,k)th  element  of  C may  be  written 


Cjk  * 1 i i n 


(A-2) 


The  covariance  matrices  associated  with  weakly  stationary 
random  sequences  are  Toeplitz  forms  with  <}>(£)  defined  by  the 
appropriate  autocovariance  or  cross  covariance  function.  The 
set  of  requirements  imposed  at  the  end  of  Section  3.2  assures 
that  the  vectors  x and  z (in  the  notation  of  the  body  of  this 
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report)  can  be  regarded  as  finite  sections  from  weakly  station- 
ary random  sequences. 

A major  step  in  the  analysis  of  frequency  domain 
methods  is  to  show  that  the  frequency  domain  representation 
of  any  Toeplitz  matrix  becomes  "asymptotically  equivalent"  to 
a diagonal  matrix  in  the  limit  as  n is  increased  without 
bound.  The  idea  of  asymptotic  equivalence  used  here  has 
previously  been  employed  in  the  analysis  of  the  limiting 
behavior  of  matrix  sequences  (Ref.  14).  This  sort  of 
special  equivalence  definition  is  valuable  because  the  element 
by  element  limit  of  the  matrices  of  interest  here  is  not , 
in  general,  diagonal.  Thus,  straight-forward  limiting 
analysis  cannot  produce  the  necessary  diagonal  estimation 
equations . 

Following  Eq . 3.2-2,  the  frequency  domain  represen- 
tation of  an  nxn  matrix  C is  defined  by 

C = FCFT  ( A-3 ) 

where  F is  the  finite  Fourier  transform  matrix  of  Eq.  3.2-3. 

The  limiting  behavior  of  C = C(n),  with  elements  defined  by 
Eqs.  A-2  and  A-3,  will  be  considered  as  n increases  without 
bound.  The  sequence  of  matrices  (C(n),  n=l,2,...}  will  be 
compared  to  a sequence  of  diagonal  matrices  (D(n),  n=l,2,...}. 
Let  the  "continued"  function  <f>c(£),  -(n-1)  <.  £ <,  n-1  be  de- 
fined by 


4»(k)  , k-0 

4>(k)  + <|>(n-k)  , k^O 


( A-4 ) 


The  discrete  power  spectral  density  of  the  matrix  C can  now 
be  defined  as 
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1 n_1 

•(«)  - = l e 

n k«0 


-2iriu)k 
n 


*c(k) 


( A-5) 


for  a)  = 0 n-1.  The  matrix  D(n)  is  defined  by 

D(n)  = n*  diag  [4>(oy , . . . ,4>(n-l)]  (A-6) 

The  asymptotic  equivalence  result  described  in  Ref.  14 
can  be  stated  as  follows 


The  matrices  C(n)  and  D(n)  are  asymptotically 
equivalent  in  the  sense  that 


Lim 
n-*-«>  n 


i E 

all  J,k 


[Cjk(n)  - Djk(n)] 


( A-7 ) 


This  is  a statement  that  the  average  squared  difference  be- 
tween C(n)  and  D(n)  vanishes  in  the  limit  as  n increases 
without  bound. 


The  asymptotic  optimality  of  the  finite,  frequency 
domain  estimation  procedure  is  established  by  examining  the 
limiting  forms  of  the  appropriate  expressions  for  rms  estima- 
tion error.  According  to  the  results  of  Section  3.2,  the 
frequency  domain  estimation  error  covariance  matrix  for  the 
optimal,  least-squares  procedure  is 


'EE 


cxx  ' cxz  czz  czx 


(A-8) 


From  Parseval's  identity,  Eq.  3.1-5,  the  rms  estimation 
error  for  the  optimal  procedure  is  given  by 
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o=  - tr( C ) 
e n v EE 


( A-9) 


The  frequency  domain  covariance  matrices  of  Eq.  (A-6)  are 
defined  by  Eq.  3.2-2.  From  Section  3.4,  the  estimation  error 
covariance  associated  with  the  finite,  frequency  domain  pro- 


cedure  is 


C'EE*  cxx  - ^zx  - cxzr'T  + gczz°T 


where  G is  the  diagonal  gain  matrix 


( A-10 ) 


G = diag 


4>  (0> 


$xz(n'1} 


*zz(n-] 


( A-ll ) 


Corresponding  to  Eq.  A-9,  the  rms  estimation  error  accompany- 
ing the  finite,  frequency  domain  procedure  is 


°e ' 3 n tr^CEE^ 


(A-12) 


Because  Eqs.  A-8  and  A-9  correspond  to  the  optimal  estimate, 
it  follows  that  oq,  2.  aQ.  The  main  result  of  this  appendix 
is  to  show  that  oe,  converges  to  the  optimal  rms  error: 


( A-13 ) 


Because  n corresponds  to  a number  of  data  points  taken  a fixed 
distance  apart,  Eq.  A-13  refers  to  a limit  as  the  survey 
track  length  is  increased  without  bound. 


The  analysis  which  establishes  Eq.  A-13  is  outlined 


as  follows: 
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First,  it  is  observed  that  Eq.  A-8  is  a 
special  case  of  Eq.  A-10,  with  gain  GQ 
given  by 

G = C C*1 
° XZ  ZZ 


( A-ll ) 


Under  the  appropriate  regularity  conditions, 
the  preceding  analysis  shows  that  the  ma- 
trices C and  C are  asymptotically  equiva- 
lent to  xz  zz 


D 

xz 


= diag[$xz(0) *xz(n-l)] 


and 


Dzz  = diae[>22(0) *zz(n-1)] 

respectively . 

• By  applying  the  theorems  of  Ref.  14  on  the 
inverse  and  product  of  asymptotically  equiva- 
lent matrices,  it  follows  that  the  gain 
matrices  G and  G0  are  asymptotically  equiva- 
lent. This  leads  to  the  conclusions  that 
CgE  and  CE£  are  asympotically  equivalent. 

• It  can  be  shown  directly  from  Eq.  A-5  that 
the  asymptotic  equivalence  of  C'  and  C--, 
implies  that  Eq.  A-10  holds. 


While  the  explicit  expressions  presented  in  this 
appendix  are  for  the  single-sensor  case  only,  the  asympto- 
tic optimality  result  of  Eq . A-13  holds  for  the  multisensor 
algorithm  as  well.  The  analysis  for  the  multisensor  case 
is  entirely  analogous  to  that  presented  here;  the  difference 
is  primarily  one  of  complexity. 


A-5 


THE  ANALYTIC  SCIENCES  CORPORATION 


APPENDIX  B 

MULTISENSOR,  DISCRETE  WIENER  FILTERING 


When  multisensor  data  is  available  t the  procedure  of 
Chapter  4 leads  to  a frequency  domain  algorithm  for  gravity 
quantity  estimation  that  incorporates  the  data  on  a sequen- 
tial, one-sensor-at-a-time  basis.  In  the  special  case  that 
all  the  data  are  available  on  "matched”  grids,  as  defined 
in  Chapter  4,  the  resulting  multisensor  algorithm  can  be  sim- 
plified into  a single,  multidimensional,  Wiener  filtering 
procedure.  The  resulting  estimation  and  estimation  error 
analysis  expressions  are  simply  vector-matrix  versions  of 
the  single-sensor  expressions,  Eqs.  3.1-3  and  3.1-4.  The  pur- 
pose of  this  apperfdix  is  to  provide  an  explicit  presentation 
of  these  expressions. 


It  is  useful  to  arrange  the  multisensor  input  data 
as  a sequence  of  vector-valued  measurements.  Let  the  measure- 
ment from  the  jib  sensor  at  the  k^h  grid  point  be  denoted 
z^,  where  1 <.  j <.  N and  1 <.  k <.  n.  Because  the  various 
sensor  data  grids  are  assumed  to  be  matched,  these  data  may 
be  arranged  in  a single,  vector-valued  sequence 


^k 


N.T 

•'Zk) 


(B-l ) 


for  k = 1 n.  The  gravity  quantity  to  be  estimated  is 

denoted,  as  before,  by  the  sequence  x^,  k = 1,  ...,  n.  The 
scalar  sequence  x^  has  the  finite  Fourier  transform  X(w), 
id  = 0 n-1,  defined  by  Eq.  3.1-1.  Similarly,  the  vector- 

valued sequence  zk  has  the  finite  Fourier  transform 
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Z(w)  = [z1(u)>...  ,ZN(0))  T]  (B-2) 

for  id  = 0 n-1,  defined  in  terms  of  the  scalar  trans- 

form of  each  of  its  elements. 

Matrix-valued  power  spectral  densities  are  required. 
The  power  spectrum  of  z^  is  the  sequence  of  N*N  matrices 
denoted  by  A (id),  u>  * 0 n-1.  The  (H,m)th  element  of 

ZZ  £ 

A ( oj ) is  the  cross  power  spectral  density  between  z (id)  and 

ZZ 

zm(w).  More  specifically,  if  <t>£m(k)  denotes  the  cross  covari- 
ance function 


wk>  ■ <B-3) 

for  1 < j,  (j+k)  < n,  then  ♦^(‘d),  defined  according  to  Eq. 

A-5  as  the  discrete  power  spectral  density  of  is  the 

(£,m)th  element  of  A . The  Nxl  cross  spectral  matrix  A (id), 

ZZ  xz 

u)  = 0,  ...,  n-1,  is  defined  in  the  analogous  element-by-element 
fashion.  As  in  the  body  of  this  report,  the  scalar  sequence 
4>  (w) , w * 0,  . . . , n-1,  denotes  the  power  spectral  density 

XX 

of  x^,  defined  according  to  Eq.  A-5. 

The  appropriate  estimation  expression  is  now  simply 

stated  as 


X (<d)  - Axz(to)  A“z  Z(io)  (B-4) 

for  id  * 0 , ...,  n-1.  The  corresponding  expression  for  the 
estimation  error  power  spectral  density  is  given  by 


4>  ( Id ) = 4>  (id)  - 

ee'  xx  ' 


Axz<“)  Azz(">  Azx(u) 


(B-5) 
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Note  that  the  matrix  inversions  indicated  in  Eq.  B-4  and  B-5 
are  of  dimension  equal  to  the  number  of  sensors,  N. 

Through  a still  more  detailed  treatment  it  is  possi- 
ble to  write  down  expressions  of  greater  generality  than  Eqs. 
B-4  and  B-5.  For  example,  these  expressions  may  be  extended 
to  the  case  where  xk  and  zfc  are  on  mismatched  grids'. 
Generalizations  of  this  type  may  be  obtained  using  the 
techniques  of  Chapter  4. 


B-3 


THE  ANALYTIC  SCIENCES  CORPORATION 


APPENDIX  C 

DERIVATION  OF  FREQUENCY  DOMAIN  STATISTICS 
USING  THE  SAMPLING  THEOREM 

The  sampling  theorem  for  finite  Fourier  transforms 
is  a useful  vehicle  for  the  calculation  of  statistics  for 
multisensor  data  on  grids  which  are  ’’mismatched"  in  the  sense 
of  Chapter  4.  Frequency  domain  autocovariances  and  cross  co- 
variances  involving  data  on  mismatched  grids  are  required  in 
the  development  of  the  estimation  procedures  of  Chapter  4. 
This  appendix  outlines  the  use  of  the  sampling  theorem  in 
deriving  the  appropriate  expressions. 

In  the  notation  of  Chapter  4,  let  u and  v denote 
data  vectors,  corresponding  to  different  sensors,  of  dimen- 
sion n and  m,  respectively.  The  information  in  u and  v is 
to  be  used  to  compute  an  estimate  of  then-vector  x.  It  is 
assumed  that  the  data  which  comprise  u and  v are  available 
on  mismatched  grids,  with  n>m.  An  estimate  for  x is  to 
be  computed  using  the  frequency  domain  version  of  the  sensor- 
sequential  estimation  expressions,  Eqs.  4.2-1  through  4.2-6. 
These  expressions  require  all  six  of  the  frequency  domain 
covariance  matrices  associated  with  the  transforms  X,  U, 
and  V i . e . , C^^  , C^y  , Cyy  , ^VX  ' ^VU  ’ * Because  x 

and  u are  of  compatible  dimension,  the  first  three  of  the 
covariance  matrices  listed  above  are  known  to  be  well- 
approximated  by  diagonal  forms  (see  Eq.  4.3-2).  The  use  of 
the  sampling  theorem  makes  it  possible  to  deduce  similar 
expressions  for  the  covariances  that  involve  the  "mismatched” 
data,  v. 
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The  elements  of  v are  to  be  viewed  as  a sampling 
of  an  expanded  data  vector  of  dimension  n,  v*.  It  is  as- 
sumed that  n = rm,  where  r is  an  integer.  The  elements  of 
v in  terms  of  those  of  v*  are  defined  by  Eq . 4.3-8.  The  re- 
sult of  the  discrete  sampling  theorem,  relating  the  elements 

£ 

of  the  m-  and  n-dimensional  transforms  V and  v , is  given 
by  Eq.  4.3-9.  This  relationship  can  be  expressed  in  vector- 
matrix  notation  as 


V = SV* 

where  S is  the  matrix 


S 


(I 


m 


m 


V 


r blocks 


(C-l) 


( C-2 ) 


The  notation  Im  indicates  the  m*m  identity  matrix. 

The  vector  v*  has  been  constructed  to  be  statis- 
tically stationary,  and  to  be  of  dimension  compatible  with 
x and  u.  It  follows  from  the  asymptotic  analysis  discussed 
in  Section  3.2  and  detailed  in  Appendix  A,  that  the  covari- 
ance matrices  Cy*U(  cy*x’  and  CV*V*  may  be  taken  to  be  diaS~ 
onal.  The  diagonal  of  each  of  these  matrices  is  the  power 
spectral  density  with  the  matching  subscript,  e.g., 

Cy*u  « diag($v*u(0) $v*u(n-l))  (C-3) 

The  required  covariance  matrices  involving  V may  now  be  ob- 
tained by  applying  Eqs.  C-l  and  C-2  to  the  corresponding  diag- 

* 

onal  matrices  involving  V . The  appropriate  relations  are 
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;vu  " sVu 


CVX  * SCV*U 


( C-4 ) 


CW  “ SCv*V*S 

The  resulting  frequency  domain  covariance  matrices  are  highly 
structured  forms  that  can  be  exploited  in  the  estimation  equa- 
tions. This  is  exemplified  by  the  diagonal-block  form 


<J>  * (0) 
v*uv  ' 


Vu(m) 


Vu(n-m) 


C 


vu 


Vu(m_1)  *vfu(2tn-l) 


$ * (n-1 ), 
v u 1 

(C-5) 


The  substitution  of  the  covariance  matrices  of  Eq.  C-4,  worked 
out  in  a form  analogous  to  Eq.  C-5,  into  Eqs.  4.2-3  through 
4.2-6  leads  to  the  estimation  equations  for  data  on  mismatched 
grids,  Eqs.  4.3-10  through  4.3-13. 


i 


* 
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